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A DESIGN HANDBOOK FOR PHASE CHANGE THERMAL 
CONTROL AND ENERGY STORAGE DEVICES 


I. INTRODUCTION 


A. Definition 

Phase change thermal control devices liave been discussed extensively in the 
literature. These articles often refer to a device of this type by different names such as 
thermal capacitor, thermal flywheel, heat of fusion device, latent heat device, and fusible 
mass device. However, all these terms refer to a component which is used to either 
thermally control a medium or store energy by utilizing a material which undergoes a 
change of phase. 

There are a number of phase change transformation classes given by Lorsch [ 1 ] as 

• Solid-liquid transformations (melting/freezing) (e.g., melting or freezing of 
ice or water) 

• Liquid-gas transformations (vaporization) (e.g., boiling of water) 

• Solid-gas transformations (sublimation) (e.g., sublimation of solid carbon 
dioxide at atmospheric pressure, “dry ice”) 

• Solid-solid transformations (e.g., transition of the rhombic torm of sulfur 
to the monoclinic form) 

• Liquid-liquid transformations (e.g., the transformation occurring when two 
liquids are immisible such as water and phenol). 

Very little energy is released or absorbed by liquid-liquid transformations, and it 
is questionable whether this is a true class of phase change. The liquid-liquid 
transformation is neglected in the remainder of this work. Liquid-gas and solid-gas 
transformation classes have large volume changes associated with them, a feature which 
requires heavy and complicated pressure vessels or special design features, such as internal 
bellows. For this reason, they also are not considered herein. Solid-liquid transformations 
are of great importance because most classes of materials undergo this type of 
transformation without exhibiting large volume changes while releasing or absorbing 
relatively large quantities of energy. The fourth transformation is also of interest because 
the energy exchange can be significant and a number of materials display this 
phenomenon in a temperature range near their melting point. Consequently, for the 
purpose of this document, phase change processes are limited to solid-liquid and 
solid-solid transformations. 



B. Applications 


There are a number of uses to which phase change devices have been applied. 
Some of these are: 

• Thermal damping of oscillatory outputs (e.g., Skylab space radiator fluid 
outlet) 

• Inhibition of thermal excursions (e.g., Lunar Roving Vehicle drive control 
electronics) 

• Maintaining constant temperature (e.g., Pegasus III coating experiment- 
thermocouple reference) 

• Energy storage (e.g., solar heating and cooling energy storage such as 
illustrated by the MIT house application) 

• Control mechanisms (e.g., Skylab 40° F vernatherm valve) 

• Temperature indication (e.g., mushroom thaw warning indicator). 

Brief descriptions of most of these applications through 1973 are given by 
Humphries [2]. Figure 1 depicts thermal responses for the first three applications. 

Typically, a phase change device is passive with no moving parts, consisting of an 
external housing and the phase change material (PCM), with or without interspersed filler 
material (Fig. 2). Usually this device is in intimate thermal contact with the medium 
which it is thermally controlling. Although shown in a simplified form, this device can be 
applied in an almost infinite number of geometrical configurations with a myriad of 
different filler configurations. 


C. Phase Change Material 

Studies have shown that to qualify as a good PCM the material should possess the 
following characteristics; 

• High heat of fusion per unit mass and volume 

• Proper melting point (or range). 

Other property requisites, but not requirements, for a good PCM are: 

• High thermal diffusivity 

• High coefficient of thermal conductivity 
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TIME 


a. THERMAL DAMPING OF OSCILLATORY OUTPUTS 



TIME 


b. INHIBITING THERMAL EXCURSIONS 



TIME 


c. MAINTAINING CONSTANT TEMPERATURE 

Figure 1. Responses of a phase change device. 
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EXTERNAL HOUSING 



• Noncorrosive 

• Low coefficient of expansion 

• High flash point 

• Good wetting characteristics 

• Minimum of void spaces (in solid-solid state) 

• Stable 

• Congruent phase change 

• Little or no supercooling 

• Relatively pure 

• Nontoxic. 

Economic factors of importance are: 

• Low cost 

• Present and future availability. 
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There are a number of classes of materials which have been investigated for use in 
phase change devices. Some of the more important are: 

• Inorganic salt hydrates (e.g., glauber’s salt Na 2 S 04 IOH 2 O) 

• Organic compounds (e.g., paraffins Cj^H 2 f,+ 2 ) 

• Eutectics of the above (e.g., 88 mole percent acetic acid + 12 mole 
percent benzoic acid) 

• Natural elements (e.g., sulphur, S) 

• Water. 

Of these materials, paraffins have been the most widely used, primarily in the space 
program. 


D. Scope of Study 

This section presents an overview of phase change devices; that is, it gives 
definitions, typical applications, and general limitations of this study. Section 11 presents 
PCM property data. Since most of the past devices have used paraffin as the PCM, 
Section II is primarily devoted to this class of materials. However, a large volume of 
nonparaffin data excerpted from other studies is given in the appendices. These data 
should allow the designer to use this document without resorting to time consuming 
literature searches for property data. Section III details the derivation and use of a 
nondimensional computer program. Section IV presents results of parametric studies 
performed using this nondimensional model. Using typical paraffin property data and 
thermal boundary conditions commensurate with earlier applications, optimum thermal 
design of the containment device is given parametrically. Examples for typical designs are 
given to illustrate use of these data. The intent of this study is to present sufficient data 
to allow the thermal designer to “short cut” time consuming analyses by either using the 
given numerical model and altering it to suit his application or, where possible, by using 
the parametric data directly. At a minimum, these data should allow the designer to get 
starting point estimates for his design. Section V outlines related topics. These topics 
include data which compare analytical and experimental results, discussing some of the 
discrepancies between the two. Section V also includes schemes for estimating convective 
effects. Other containment cases not covered by the examples used in Section IV are 
also given in Section V. 

The reader should keep in mind that the treatment of phase change devices given 
in this document is limited to thermal aspects. As such, other topics of importance such 
as structural and material designs and selection are not addressed. Also, specialized 
subject matter such as the kinetics of the phase change process and chemical 
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plienomenon (e.g., supercooling, incongruent melting, and dendrite formation) are not 
discussed. Although convection analysis is discussed, details of convection stimulii such as 
buoyancy and surface tension forces are ignored. The reader is referred to survey reports 
by Grodzka [3,4] and a text by Chalmers [5] for information in these areas. Due to the 
large number of possible design configurations, this report addresses only typical aird 
simplified schemes. 


II. PROPERTY DATA 


In several space-related applications of phase change devices, a paraffin has served 
as the PCM. Major emphasis in this section is given to those properties which are 
typically needed in a thermal design of selected paraffins. These include the melting 
temperature, transition temperature, heat of fusion, heat of transition, specific heat, 
density, thermal conductivity, viscosity, and coefficient of volume expansion. 

Since hydrocarbon properties are dependent on purity, differences in reported 
property values may be directly related to differences in purity, a quantity often not 
specified. 

Although this section is devoted to properties of paraffins, some of the properties 
of certain nonparaffins are presented in Appendix A. These have been extracted from 
several different references, and, in some cases, some paraffins and paraffin mixtures are 
included there also. 

It is understood that the term paraffin generally denotes any of the saturated 
aliphatic hydrocarbons of the methane series Property data herein are 

principally limited, however, to those members of the family lying between n-Undecane 
(C 11 H 24 ) and n-Triacontane (C 30 H 62 ) because the corresponding range on melting 
temperature has included that of previously considered applications. 


A. Sources, Availability, and Estimated Cost of Selected Paraffins 

Table 1, which is based on a reported survey [6], shows a limited selection of 
companies and laboratories capable of supplying several selected paraffins in quantities 
adequate to be used in phase change devices. Also, the availability of Sunoco P-116 
paraffin wax from the Sun Oil Co., Ci 4 -C ,6 paraffin from Conoco, and a C, 5 -Ci 6 
paraffin mixture from Enjay is noted in Reference 1. 

Information relative to type, grade, and physical property data availability from 
the suppliers listed in Table 1 is outlined in Table 2 based on a review [6] of their 
products. 
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TABLE 1. POTENTIAL SUPPLIERS OF SELECTED PARAFFINS^ 


Company 

Address 

American Petroleum Institute 
Standard Reference Materials 

Carnegie-Mellon University 
Schenley Park 

Pittsburgh, Pennsylvania 15213 

Aldrich Chemical Co., Inc. 

940 West Saint Paul Avenue 
Milwaukee, Wisconsin 53233 

Analabs Inc. 

80 Republic Drive 

North Haven, Connecticut 06473 

Chemical Sample Co. 

4692 Kenny Road 
Columbus, Ohio 43221 

Eastman (Kodak) Organic 
Chemicals 

2400 Mt. Read Boulevard 
Rochester, New York 14650 

Gallard-Schlesinger Chemical 
Manufacturing Corp. 

584 Mineola Avenue 

Carle Place, New York 11514 

The Humphrey Chemical Co. 

Devine Street 

North Haven, Connecticut 06473 

Lachat Chemicals Inc. 

20200 Ashland Avenue 
Chicago Heights, Illinois 604 1 1 

Phillips Petroleum Co. 

Bartlesville, Oklahoma 74004 

Polysciences Inc. 

Paul Valley Industrial Park 
Warrington, Pennsylvania 18976 


a. From Table 9 of Reference 6. 


The cost of paraffins is influenced significantly by purity, and it also varies 
considerably with source. The ultrahigh purity standard reference materials available from 
the American Petroleum Institute (API) cost $95 for a 5 ml unit [6]. Some price 
estimates for two grades of several paraffins are given in Table 3. All but the last three 
entries in Table 3 correspond to prices for relatively small quantities of reasonably high 
purity paraffins, probably made available for laboratory rather than large-scale 
commercial purposes. The last three entries denote some paraffin mixtures which can be 
obtained in large quantities. 
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TABLE 2. CAPABILITIES OF POTENTIAL SUPPLIERS OF SELECTED PARAFFINS^ 
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From Table 10 of Reference 6. 




TABLE 3 . SOME APPROXIMATE PRICE DATA FOR SELECTED PARAFFINS 


Material 

Formula 

Purity 

(percent) 

Approximate^ 
Unit Cost 
$/kg ($/lb) 

Quantity 
kg (lb) 

Reference 

n-Nonane 

C9H20 

99 

26.46 (12.00) 

2.72 (6) 

6 

n-Nonane 

C9 0 

95 

14.29 (6.48) 

2.72 (6) 

6 

n-Decane 

C,oH22 

99 

16.26 (7.38) 

2.77 (6.1) 

6 

n-Decane 

Cj 0^22 

95 

12.87 (5.84) 

2.77 (6.1) 

6 

n-Undecane 

Cj 1 H24 

99 

23.47 (10.65) 

2.81 (6.2) 

6 

n-Undecane 

C,iH24 

95 

14.21 (6.44) 

2.81 (6.2) 

6 

n-Dodecane 

Ci2H,6 

99 

16.45 (7.46) 

2.86 (6.3) 

6 

n-Dodecane 

Cj 2H26 

95 

12.65 (5.74) 

2.81 (6.2) 

6 

n-Tridecane 

Cl 3H28 

99 

34.99 (15.87) 

2.86 (6.3) 

6 

n-Tridecane 

Cl 3H2 8 

95 

22.29 (10.11) 

2.86 (6.3) 

6 

n-Tetradecane 

Cl 4H30 

99 

18.95 (8.59) 

2.90 (6.4) 

6 

n-Tetradecane 

Cl 4H30 

95 

15.62 (7.09) 

2.90 (6.4) 

6 

n-Pentadecane 

C,sH32 

99 

41.34 (18.75) 

2.90 (6.4) 

6 

n-Pentadecane 

C,sH32 

95 

30.13 (13.67) 

0.68 (1.5) 

6 

n-Hexadecane 

Cl 6^34 

99 

35.27 (16.00) 

0.68 (1.5) 

6 

n-Hexadecane 

Cl 6^34 

95 

20.08 (9.11) 

2.90 (6.4) 

6 

n-Heptadecane 

Cl 7H36 

99 

37.89 (17.19) 

2.90 (6.4) 

6 

n-Heptadecane 

Cl 7^36 

95 

30.13 (13.67) 

0.68 (1.5) 

6 

n-Octadecane 

Cl 8^3 8 

99 

18.65 (8.46) 

2.95 (6.5) 

6 

n-Octadecane 

C18H38 

95 

30.13 (13.67) 

0.68 (1.5) 

6 

n-Nonodecane 

Cl 9 H4 0 

99 

40.70 (18.46) 

2.95 (6.5) 

6 

n-Nonodecane 

Cl 9H4 0 

95 

27.93 (12.67) 

0.68 (1.5) 

6 

n-Eicosane 

C2 0 H4 2 

99 

26.05 (11.82) 

2.99 (6.6) 

6 

n-Eicosane 

C20H42 

95 

41.15 (18.67) 

0.68 (1.5) 

6 

P -116 Paraffin 


- 

0.14 (0.065) 

Not specified 1 

Ci5“Ci 6 Paraffin mixture 

- 

0.11 (0.05) 

4536 (10,000) 1 

Ci4-Ci 6 Paraffin 

1 

- 

0.11 (0.05) 

Not specified 1 

1 


a. Prices shown are approximate based on 1974 publications. All but the last three 
entries shown in the table correspond to prices quoted for the small quantities 
indicated which, more than likely, are for laboratory applications and are not 
totally indicative for large commercial orders. Also, it should be noted that 
there is also variation in prices with supplier. 
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B. Melting and Transition Temperatures 


Table 4 contains published [7,8] values of the melting and transition 
temperatures for the listed paraffins. 

For commercial grades of paraffins, some variations in the phase change 
temperatures are to be expected. Impurities even tend to cause some difference in the 
melting and freezing temperatures. Some reported [ 6 ] experimental results which show 
the effect of purity on the fusion temperature for several grades of certain paraffins 
supplied by two manufacturers are given in the lower two segments of Table 5 together 
with the corresponding fusion temperature for the high purity paraffins available from 
API. Similar results and API data for solid-solid phase transition temperatures are given in 
Table 6 . 

For the seven paraffins listed in Tables 5 and 6 , correlations of phase change 
temperatures corresponding to API data with the number of carbon atoms in the chain 
are shown in Figure 3. 


C. Heat of Fusion and Heat of Transition 

The energy values associated with solid-liquid and solid-solid transitions as 
published by API [7] or either by Timmermans [ 8 ] are also listed in Table 4. The 
effects of purity on these energy values based on comparisons between measurements [ 6 ] 
on two commercial grades (lower two segments of Tables 7 and 8 ) and values reported 
by the API for high purity paraffins are listed in Tables 7 and 8 , respectively, for several 
paraffins. 

For the seven paraffins listed in Tables 7 and 8 , correlations of the heat of fusion 
and the heat of transition with the number of carbon atoms in the chain are shown in 
Figures 4 and 5, respectively, for the data available from API. 

It is interesting to note also the effect of the mixing of two different paraffins on 
the heats of fusion and transition. As an example, the results reported [ 6 ] for a 
hexadecane-octadecane binary system are shown in Figure 6 . 


D. Specific Heat 

The specific heat is temperature dependent in each phase with an increase in 
specific heat occurring with an increase in temperature. There is a change in specific heat 
between two phases. Data for all the paraffins listed in Table 4 have not been obtained. 
For some of the listed paraffins, specific heat data over a range of temperature are given 
in Reference 8 . These have been plotted in Figures 7, 8 , and 9. There are obvious gaps in 
the data in the vicinity of phase transitions. Data given in Reference 9 for C, 4 H 3 q, 
^ 16 ^ 34 , C, 8 H 38 , and C 20 H 42 are shown in Figures 10 through 13, respectively. 
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TABLE 5. FUSION TEMPERATURE FOR SELECTED PARAFFINS^ 



* 
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Based on Tables 1 and 2 of Reference 6. 

* Supplier A 

* Supplier B 

* Supplier C 




TABLE 6. TRANSITION TEMPERATURES FOR CERTAIN PARAFFINS^ 



XU 

H < « U 


^ CXi cx 

^ ^ & 9^ 

^ D 

^ CO c/3 c/5 
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PHASE TRANSITION TEMPERATURE, 



Figure 3. 


Correlation of phase change temperature with number of carbon atoms 
in the chain (from Fig. 1 of Reference 6). 
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PHASE TRANSITION TEMPERATURE, 



TABLE 7. COMPARISON OF HEAT OF FUSION FOR SEVERAL 
GRADES OF SELECTED PARAFFINS 


Material 

Heating 

(Solid Liquid) 
kJ/kg (Btu^b) 

Cooling 

(Liquid Solid) 
kJ/kg (Btu/lb) 

Purity 
(Mole %) 
or Grade 

n-Undecane 

142 ( 61.2) 

142 ( 61.2) 

99.96 

± 0.03 


n-Dodecane 

216 ( 93.2) 

216 ( 93.2) 

99.969 ± 0.025 


n-Tridecane 

154 ( 66.5) 

154 ( 66.5) 

99.91 

± 0.06 


n-Hexadecane 

235 (101.5) 

235 (101.5) 

99.90 

± 0.06 


n-Octadecane 

244(105.0) 

244 (105.0) 

99.90 

± 0.08 

1 

n-Nonadecane 

187 ( 80.6) 

187 ( 80.6) 

99.90 

± 0.08 


n-Eicosane 

248 (107.0) 

248(107.0) 

99.90 

± 0.08 

J 


n-Undecane 

140 ( 60.4) 

143 ( 61.7) 

99 ^ 



n-Dodecane 

218 ( 94.0) 

220 ( 94.9) 

99 



n-Tridecane 

145 ( 62.5) 

159 ( 68.5) 

99 



n-Hexadecane 

237 (102.0) 

235 (101.5) 

99 


n-Octadecane 

248 (107.0) 

245 (105.6) 

99 



n-Nonadecane 

180 ( 77.5) 

180 ( 77.5) 

99 



n-Eicosane 

249 (107.5) 

245 (105.5) 

99 ^ 



n-Undecane 

142 ( 61.4) 

141 ( 61.0) 

99.75 



n-Dodecane 

218 ( 94.0) 

212( 91.5) 

99 



n-Tridecane 

154 ( 66.2) 

150( 64.8) 

99 



n-Hexadecane 

239 (103.0) 

235 (101.5) 

99 



n-Hexadecane 

244(105.0) 

237 (102.0) 

95 



n-Octadecane 

243 (104.6) 

242(104.5) 

95 



n-Nonadecane 

182 ( 78.3) 

179 (77.2) 

95 



n-Eicosane 

255 (110.0) 

243 (104.8) 

95 



n-Eicosane 

241 (104.0) 

234(101.0) 

90 

> 



* Supplier A 
Supplier B 
*** Supplier C 


15 



TABLE 8. COMPARISON OF HEAT OF TRANSITION FOR SEVERAL GRADES OF 

SELECTED PARAFFINS 


Material 

Heating 
(Solid Solid) 
kJ/kg (Btu/lb) 

Cooling 
(Solid ->■ Solid) 
kJ/kg (Btu/lb) 

n-Undecane 

n-Tridecane 

n-Nonadecane 

439 (18.9) 
416 (17.9) 
514 (22.2) 

439 (18.9) 
416(17.9) 
514(22.2) 

n-Undecane 

n-Tridecane 

n-Nonadecane 

409(17.6) 
433 (18.7) 
513 (22.1) 

405 (17.4) 
470 (20.3) 
516(22.3) 

n-Undecane 

n-Tridecane 

n-Nonadecane 

403 (17.4) 
403 (17.4) 
495 (21.4) 

422(18.2) 
406(17.5) 
499 (21.5) 


Supplier A 
Supplier B 
Supplier C 


Purity 
(Mole %) 
or Grale 


99.96 

99.91 

99.90 













HEATS OF FUSION AND TRAN^TIOFMcJ/kS 



10 15 20 

CARBON ATOMS IN CHAIN 


Figure 5. Correlation of heat of transition with number of carbon atoms in the 
chain [from Fig. 5(b) of Reference 6]. 



^16*^34 WEIGHT PERCENT OCTADECANE 

HEXADECANE 


CtS^Sfl 

OCTADECANE 


Figure 6. Example illustrating effect of composition of heats of phase change for 
hexadecane-octadecane binary system (from Fig. 29 of Reference 6). 
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SPECIFIC HEAT, kJ/t<g-®C 


TEMPERATURE. “F 





Figure 8. Specific heat versus temperature for C14H30, C15H32) 3nd C15H34 
(based on data given in Reference 8). 
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SPECIFSC HEAT, Btu/Ib-®F 





TEMPERATURE. ®F 



30 40 50 60 70 80 90 100 110 

TEMPERATURE, <*C 


Figure 9. Specific heat versus temperature for C 2 sHjj 
(based on data given in Reference 8 ). 


Measurements of specific heat for several paraffins using samples from two 
manufacturers are reported by Reference 6 . These are tabulated in Table 9 together with 
values for the same materials from API. 

Comparative examination of data shown in Figures 7 through 13 and listed in 
Table 9 indicates that the specific heat of paraffins is sensitive to impurities. 


E. Density 

The density of the liquid phase decreases with temperature and data for most of 
the considered paraffins are shown in Figure 14. The values agree with those reported in 
Reference 6 . For Ci 4 H 3 o> C 14 H 34 , CigH 3 g, and C 20 H 42 , density data are also plotted 
in Figures 15 through 18, respectively, which are from Reference 9 and which include 
some solid phase density measurements. For these paraffins, density changes across the 
solid-liquid phase transition ranged from 5 to 8 percent. 


F. Thermal Conductivity 

There is a limited amount of thermal conductivity data available for paraffins. 
Most reported values are only for the liquid phase. For selected paraffins Table 10 lists a 
limited amount of liquid thermal conductivity data, which were reported in Reference 6 
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Figure 10 . Specific heat of tetradecane (C14H30) (from Fig. 8 of Reference 9 ) 
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trence 9). 



TEMPERATURE 



Figure 13. Specific heat of eicosane (CjoH 4 j) (from Fig. 1 1 of Reference 9). 




TABLE 9. SPECIFIC HEAT DATA FOR CERTAIN PARAFFINS^ 
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TEMPERATURE, ®P 



TEMPERATURE. ®C 


Figure 15. Density ot tetradecane <Ci 4 H 3 o) (from Fig. 4 of Reference 9) 


TEMPERATURE, ®F 



TEMPERATURE, °C 


Figure 16. Density of hexadecane (C 16 H 34 ) (from Fig. 5 of Reference 9) 






TABLE 10. THERMAL CONDUCTIVITY OF SELECTED PARAFFINS (EVALUATED 

AT PHASE CHANGE TEMPERATURE)^ 


Material 

W/m-°C 

Btu/h-ft-°F 

n-Undecane 

0.1496 

0.0865 

n-Dodecane 

0.1488 

0.0860 

n-Tridecane 

0.1496 

0.0865 

n-Hexadecane 

0.1505 

0.0870 

n-Octadecane 

0.1505 

0.0870 

n-Nonadecane 

- 

— 

n-Eicosane 

0.1505 

0.0870 


a. From Table 8 of Reference 6 . 


as being available from the API. Thermal conductivity data which are presented in 
Reference 9 for the liquid phase of C 14 H 30 , Cj 6 H 34 , Ci gHag, and C2oH42> ^re shown 
in Figures 19 through 23. Figure 23 shows a correlation which was extrapolated to 
obtain thermal conductivities for C 2 oH 4 2 - 

The lack of thermal conductivity data for the solid phase poses some analytical 
uncertainties when attempting to analyze the phase change process. This concern and 
related considerations are treated in more detail in Section V. 


G. Viscosity 

Absolute viscosity data for liquid paraffins C 11 H 24 through C 20 H 42 . which are 
given in Reference 7, are shown plotted as a function of temperature in Figure 24. 


H. Surface Tension 

Values of the surface tension for liquid paraffins C, ,H24 through C 20 H 42 , which 
are tabulated in Reference 7, are plotted versus temperature in Figure 25. The data varies 
linearly with temperature with an increase of temperature resulting in a decrease of 
surface tension. 


I. Coefficient of Expansion 

Single values of the coefficient of expansion for four of the considered paraffins 
are shown in Table 1 1 which were reported in Reference 6 as being from API. 
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Figure 20. Thermal conductivity of hexadecane (C , 41134 ) 
(from Fig. 13 of Reference 9 ). 
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Figure 22. Thermal conductivity of eicosane (C 20 H 42 ) extrapolated from data 
for lower carbon number paraffins (from Fig. 15 of Reference 9 ). 
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Figure 23. Correlation of thermal conductivity data for normal 
paraffins (CnH2n+2^ Reference 9 ). 
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TABLE 11. VOLUMETRIC COEFFICIENT OF EXPANSION (60°F/15.6^C)^ 


Material 

l/K 

l/^F 

n-Undecane 

10.1 X 10'“* 

5.6 X 10*'’ 

n-Dodecane 

9.9 X 10“* 

5.5 X lO*'’ 

n-Tridecane 

9.4 X lO’'* 

5.2 X lO"* 

n-Eicosane 

8.5 X lO-'* 

4.7 X 10*'’ 


a. From Table 18 of Reference 6 . 


ill. TWO-DIMENSIONAL ANALYSIS OF A PHASE 
CHANGE DEVICE 


A. introduction 

The analysis of two-dimensional conductive heat transfer within a phase change 
device is outlined in this section. Since convection may be present in certain applications, 
the design of a phase change device based solely on conduction will in most cases be 
conservative (i.e., heat transfer rate to and from material will be lower than in a case 
where convection is present). Additional attention is given in Section V to the effect of 
convection and how it can be included in the analysis. 

The two-dimensional transient heat conduction equation for an isotropic, 
homogeneous medium in the absence of sources, sinks, or phase change is 



4 

In general, solutions of equation (1) for typical applications can best be obtained 
via numerical techniques, other techniques were outlined by Muehlbauer and Sunderland 
[10]. For numerical solutions of equation (1), the partial derivatives are approximated by 
finite difference analogs. A forward difference analog for the time derivatives can be 
used. With regard to the latter, there are two possibilities for a transient problem [11]. If 
all temperatures appearing in the spatial analogs are evaluated at the old time level, the 
formulation is explicit, and the resulting equation for a node involves only one unknown 
at the new time level. The solution for a system of nodes by this method is fairly simple 
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and straightforward. If all temperatures appearing in the spatial analogs are evaluated at 
the new time level, the formulation is implicit. The resulting equation for a node involves 
five unknown temperatures at the new time level. Solution for a system of nodes involves 
solving a set of simultaneous equations equal to the number of nodes. 

Figure 26 represents a general node (I, J), its surrounding nodes, and the 
notational scheme for interconnecting thermal conductances. 



Figure 26. General designations for nodes and connecting conductances. 


Physical arguments can be used to formulate the difference analogs previously 
discussed. The summation of heat transfer to node (I, J) from all surrounding nodes 
should equal the product of the nodal capacitance and the difference in its old and new 
temperature. With reference to Figure 26, the explicit formulation based on this physical 
approach is 

GH(I,J)[T(I-1,J)-T(I,J)] +GH(I+1,J)[T(I+1,J)-T(I,J)] + GV(I,J)[T(I,J-1) 

-T(I,J)] +GV(I,J+1)[T(I,J+1)-T(I,J)] = — (2) 
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The implicit formulation is 


GH(1,J) [f (I-U)-f (U)] +GH(I+1,J)[T(I+1,J)-T(U)] +GV(U) [T(U-D 
- f (I,J)1 + GV(I,J+1 )[f (U+1 ) - T (I,J)] = £(I.J)[T (I,J) - T(I,J)1 ^ 


The tildes denote temperatures at the new time level. 


Numerical solution of the explicit formulation given by equation (2) is subject to 
a restriction on the allowable time step given by 


At < 


C(U) 

GH(I,J) + GH(I+1,J) + GV(I,J) + GV(I,J+1) 


(4) 


For a nodal arrangement, the node having the smallest time step given by equation (4) 
governs the entire solution. 

The implicit formulation given by equation (3) is not restricted by a stability 
requirement, but the solution is not as simple as one using the explicit formula. In 
general, it requires a matrix inversion routine and may be impractical for a large number 
of nodes. 

A technique for handling the phase change when using the explicit formulation 
involves keeping a record of the stored energy at each PCM node. The stored energy at 
any time is compared to that associated with the initiation and termination of the phase 
change. When the accumulated energy at a node lies between these two values, the 
temperature prediction is overridden and the temperature is forced to be equal to the 
phase change temperature. 

The technique for handling the phase change when using the implicit formulation 
involves comparison of the predicted temperature with the phase change temperature. 
When the predicted temperature first exceeds the latter value, the predicted value is 
overridden and the nodal temperature is specified to be the phase change temperature. 
The energy associated with the excess of the predicted value over the phase change 
temperature is calculated and allowed to accumulate. Once this accumulated energy 
reaches a value sufficient to accomplish the phase change, the predicted temperature is 
allowed to rise in accordance with the governing conduction equations. 
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B. Physical Model 


The physical model selected for study corresponds to a phase change device 
appUcation which is depicted in Figure 27. Basically, the phase change device consists of 
a metallic housing in the shape of a parallelepiped which is filled with PCM and straight 
metalUc fins which are included to improve the heat transfer mechanism between the 
housing and the PCM. 
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Figure 27. Phase change device application selected for study. 


When the phase change device is operating as an energy sink (i.e., it is storing 
energy), either a hot fluid or a hot surface is considered to be in contact with one 
surface of the housing resulting in heat transfer into the unit. Hereafter, this surface 
which is in contact with the hot medium is referred to as the base. The PCM is used to 
absorb the energy primarily via the phase change process. The stored energy may be later 
rejected by means of exposing the base or the cover or both to a cooler medium, which 
in some applications may be the same medium from which the energy was absorbed 
earlier. 


A typical single cell of the device is shown in Figure 28. It is assumed that planes 
parallel to the fins and passing through their centers represent planes of symmetry as do 
planes parallel to the fins and located midway between two adjacent fins. Consequently, 
the analysis is reduced to that of a half-cell represented by the schematic given in Figure 
29. 
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PLANES OF SYMMETRY 



Figure 28. Typical single cell of phase change device. 


The following thermal analysis pertains to the symmetrical half-cell shown in 
Figure 29 with the base exposed to a constant heat flux density. The resulting computer 
program is also modified for the case of a constant base temperature. 

The limitations of the analysis are given by the following assumptions: 

• Conduction is the only heat transfer mechanism within the 
two-dimensional half-cell; the overall phase change device is considered to 
be composed of many such cells. (The effects of convection are discussed 
in Section V). 

0 problem is two-dimensional; the cell, therefore, is sufficiently long 

that end effects are small. 

• The cover is perfectly insulated. 
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Figure 29. Symmetrical half-cell used in analysis. 

• The physical properties, density, specific heat, and thermal conductivity of 
the metallic fm, cell base, and cell cover are considered to be identical and 
to have constant values. Likewise, the density, specific heat, and thermal 
conductivity of the soUd and liquid phases of the PCM are considered to 
be identical and to have constant values. (The model may be easily 
modified to input temperature variant properties.) 

To cast the problem in a nondimensional formulation, appropriate dimensionless 
parameters are sought. As a means of obtaining the pertinent dimensionless parameters, 
differential equations for the temperature in the PCM, the fin, the base, and the cover are 
presented and nondimensionalized. The differential analysis is followed by a numerical 
analysis which is used for the solutions, the results of which are presented in Section IV. 
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C. Differential Formulation 


Since the principal objective of the differential formulation is to find the 
appropriate dimensionless parameters for the physical model previously described, the 
reader interested only in using the results can go directly to the next section on the 
numerical analysis. 

The differential equations pertinent to each region of the model are presented 

below. 


1. PCM Region. The energy differential equation is the usual two-dimensional 
conduction equation 


a^T _ 1 9T 

ax^ ay^ a at 


and the two-dimensional interfacial boundary condition according to Rathjen and Jiii 
[12] is 


k 



as 

phf- 


y = 6 


( 6 ) 


where 5(x,t) is the interfacial location or fusion front. Choosing the width W of the PCM 
region shown in the symmetrical half-cell in Figure 29 as the characteristic x-dimension, 
the PCM height H as the characteristic y-dimension, the PCM fusion temperature as the 

characteristic temperature, and denoting a characteristic time by t^. result in the following 
dimensionless quantities; 


x' = x/W 


(7) 


y' = y/H 


( 8 ) 


d = T/Tf 


(9) 


t' = t/t 


c 


( 10 ) 
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In terms of these, equation (5) can be written as 


^ ^ 9y'^ ate 


and equation (6) becomes 


1 

[ (-S 

hfH= /; 

Wl, 

[ \w/ \ax7_ 

oc Cp Tf V 


( 11 ) 


( 12 ) 


2. Fin. Treating the temperature distribution in the fin as one-dimensional, i.e., 
Tm = Tn^(y.t). with conductive heat transfer from the fin to the PCM, which is analogous 

to the treatment of convective heat transfer in the conventional one-dimensional fin 
problem, results in 


l‘„\S2A»Vo.y at 

where the coordinate system is given in Figure 30. The second term, on the left in 
equation (13) results from considering conduction in the PCM. The temperature in the 
fin at X = 0, however, is the same as that in the PCM at x = 0. 



Figure 30. Coordinate system and dimensions for fin. 
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Nondimensionalization of equation (13) using equations (7) through (10) results in 


^ \S2W/\ax7o y (^atJ\aj„j 9t' 


(14) 


3. Cell Base. Treating the cell base as a lumped system, i.e., Tjj = Tjj(t), and 

assuming the base and fin properties to be the same result in the following energy 
equation: 







dTb 

dt 


(15) 


where it has been assumed that the base is exposed to a constant heat flux density, q”. 
The coordinate system and physical dimensions for the cell base are given in Figure 31. 
The subscripts x,0 on the third term on the left side of equation (15) denote evaluation 
at 0.5 S 2 < X < W and y = 0. In terms of the previously defined dimensionless variables, 
equation (15) becomes 



(16) 

In equations (15) and (16) the overbar on the third term On the left side denotes a value 
spatially averaged over the x-length W. 



Figure 3 1 . Coordinate system and dimensions for base. 


45 


4. Cell Cover . Since the treatment of the cell cover is analogous to that of the 
cell base except for the assumption that the cover is insulated, no new dimensionless 
parameters should appear; therefore, the appropriate equation is omitted for brevity. 

To this point the dimensionless equations contain a undefined reference time, t^,. 
By selecting this characteristic time to be given by 


H" 



the ratio /cct^ becomes 1.0 in equations (11), (12), (14), and (16). Examination of the 

reduced set of equations reveals seven salient dimensionless ratios for this class of 
problem. Summarizing these ratios by category, with pure constants omitted, yields 


Geometric: 


R| - H/W 
R2 = H/Sj 
RS = W/S2 


Material Properties: R^ = 


■'m 

IT 


Rc = 


Ra = 


ocja 

hf 


Thermal Loading: Ry = 3 — 

k Tf 


These ratios developed with the differential equation formulation of the problem are also 
appropriate for the finite difference formulation. Hence, the finite difference equations 
are cast in dimensionless form with the seven aforementioned dimensionless parameters 
appearing explicitly and serving as focal points for a parametric numerical study. 
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D. Finite Difference Formulation 


The finite difference formulation begins with a subdivision of the region of 
interest into a network of nodal points. For the two-dimensional, symmetrical half-cell 
shown in Figure 29, the nodal network and nodal point designation scheme are shown in 
Figure 32. The approach provides flexibility for increasing or decreasing the number of 
PCM nodes in accordance with the dictates of size and shape, there being M by N PCM 
nodes. Also note that all base nodes correspond to J = 1 and all fin nodes correspond to 
1=1. The number of base or fin nodes is, of course, directly related to the PCM 
subdivision. Figure 26 depicts a general node (1,J), its four surrounding nodes, and the 
notational scheme for interconnecting thermal conductances. 



Figure 32. Nodal grid designation. 


Physical arguments can be used to formulate the temperature finite difference 
equations. With reference to Figure 26 the net heat transfer to node (I,J) from its 
surrounding nodes (during time At) should equal the product of the nodal thermal 
capacitance and the temperature difference occurring during the time increment. The 
explicit formulation based on this approach in nondimensional form is 
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GH'(U) [0(1-1, J)-0(U)] + GH'(l+l,J)t0(I+I,J)-0(U)l + GV'(I,J)[0(I,J-1) 


-0(I,J)J + GV'(I.J+1) [0(U-H)-0(I,J)] = 


CUJ)[0(U)-0(I.J)] 

M* At' 


(18) 


The implicit formulation is 


GH'(U) (0(1-1, J)-0(U)] + GH'(I-H,J)[0(I-H,J)-0(I,J)1 + GV'(U)[0(1,J-1) 


-0(1,J)] GV'(I,J+1) [0(I,J+1)-0(U)1 = 


C'(I,J) (0(I,J) - 0(I,J)1 
R,* M* At' 


(19) 


In the nondimensionalization of all thermal conductances, the nondimensional value is 
obtained by dividing the respective conductance by kB which is the conductance between 
two general PCM nodes. Similarly, nondimensional thermal capacitance values are 
obtained by dividing the respective dimensional value by pS^BCp which is the thermal 

capacitance of a general PCM node. It should be noted that the nondimensional 
conductances and capacitances that are entirely within the PCM have the numerical value 
of 1.0. All other dimensionless conductances and capacitances have values other than 
unity. As an example, the dimensionless horizontal conductance between the fin and the 
lower left hand corner PCM node is 


GH'(2,2) 


2 

2R3R4 


( 20 ) 


A dimensionless base nodal capacitance, for example, is 


C'(2,l) 


M RiR4 

Rj R5 


( 21 ) 


1. Computer Programs . Numerical solutions of the problem described herein via 
computers can be affected by use of either the explicit formulation given by equation 
(18) or the implicit formulation given by equation (19). The two general approaches 
together with the inherent advantages and disadvantages related to each are discussed in 
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# « 

the general heat transfer literature; see, for example, Chapman [11] or Ozisik [13]. It is 
well established that solution of the explicitly formulated set of equations is subject to a 
stability restriction which dictates the maximum allowable time increment. This 
restriction is given by equation (4). For a small grid size, this requirement leads to 
excessive computer time. While solutions of the implicitly formulated set of equations is 
not restricted by a stability requirement, the method involves the solution of a set of 
simultaneous equations equal in number to the number of nodes. For a small grid size, 
large computer storage is needed. 

The problem under consideration involves thin metallic fins; most attention, 
consequently, has been given to the implicit approach. The stability restriction dictated 
by a node located in the fin resulted in extremely small time increments for the explicit 
approach. 

A computer program based on the implicit approach, equation (19), for the 
problem described in this section is given in Table 12. A modification of the program to 
accommodate a sudden change in the base temperature rather than imposition of a 
constant heat flux density is given in Table 13. In both programs, a band algorithm 
technique given by Funderlic [14] is used to solve the equations. A dimensional 
computer program based on the explicit approach is presented and briefly discussed in 
Appendix D. 

2. Instructions for Using Computer Programs. To use the programs, the following 
steps are suggested: 

a. Establish the grid size (see Fig. 32) by specifying M and N. M represents the 
number of vertical columns and N represents the number of horizontal rows within the 
PCM. The ratio of N to M must equal Rj. The accuracy of the transient results may be 
affected by the choice of M and N, particularly in the early part of the transient. The 
user should conduct a sensitivity study by altering M and N along with DT (Step f) in a 
systematic manner until the results appear to be unaltered to within an acceptable level 
by further changes. The influence of the choice of M and N is discussed and illustrated in 
Section IV along with discussion of the parametric study. 

b. Determine appropriate values for the dimensionless ratios Rj through R 7 as 
defined in the analysis. Note that R 7 is denoted as Rg in the program. Also, for the 
constant temperature case R 7 is not needed, but a temperature difference (DELT) is 
needed. This difference is defined by 

DELT = 1 

Tf 

c. The dimension card must be set correctly. The dimensions depend on M and 
N as indicated below. 


A(LL,L1) 

T1(M1,N2) 


49 



T2(LL) 

GH(M2,N2) 

GV(M1,N2) 

C(M1,N2) 

F(LL) 

TEXC(LL) 

D(LL) 

MELT(LL) 


where 


LL = (N2) X (Ml) 

LI = (2) X (Ml + 1) 

Ml = M + 1 
M2 = M + 2 
N2 = N + 2 

d. Input required parameters in accordance with Formats 998 and 999. The 
required inputs are shown by the Read statements. NR is reserved for a reference run 
number. 


e. The program is written with the initial PCM temperature at Tp If it is desired 
to run for other values of initial temperature, set TIN by 


TIN = Tj^itial/Tf 

f. Select a nondimensional time set by specifying DT. It is suggested that the 
magnitude of this time increment be explored for its influence on the accuracy of the 
results. Note the discussion given in Step a also. 

g. Set ITT. This quantity controls the printout of results. The program is written 
to yield output at the end of the first time increment and then on intervals given by 


(ITT - 1 ) X (DT) 
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The output includes 


TIME 

TRISE 

UCOEF 


- Nondimensional time measured from start of melting 
-- (Tg^gg - Tj-)/Tj- (use absolute temperatures) 

-- Nondimensional overall heat transfer coefficient defined by 


UCOEF 




Tease -Tf 


FMELT — Fraction of PCM melted at given time (or solidified for freeze ease) 

SQF - Energy stored in fin 
SQB - Energy stored in base 

SQT -- Energy stored in top 

SQS Energy stored in solid 

SQL -- Energy stored in lujiiid 
QLAT Energy used to change phase 

QTOT — Summation of SQF, SQB, SQT, SQS, SQL and QTAF 

EIN - Total energy transferred into base which is dependent on q” and time 

RATIO - QTOT/EIN check on energy balance. 

All the energy values which were listed are nondimensional and evaluated with respect to 
initial values of zero. The energy value used in the nondimensionalization is the latent 
energy for a single PCM node. 

The programs given in Tables 12 and 13 can be altered slightly to accommodate 
freezing for the case of a constant heat flux density withdrawal at the base and a step 
reduction in the base temperature, respectively. The changes needed to accomplish this 
are; 


a. In the parameters at the beginning, replace QLATM by -QLATM. 

b. In the third, second, and first lines above statement 150, change the 
statements to 
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TABLE 1 2. COMPUTER PROGRAM FOR BASE EXPOSED TO A CONSTANT 

HEAT FLUX DENSITY 


c constant heat rate program 

OIhENSION A(1E6«7)« TlO«Aa>« T2(126)i 3H(4«A2)« QVO 
1«42)«C(3M2)4 P(126)# TEXC(126M DU26I# iELT|126) 
IMPLICIT REAL»8( A-H#0-Z» 

Rl • H/W ; R2 ■ H/Sl 1 R3 ■ W/S2 J R4 - </iO< / «««*««« 

R5 ■ (K/DEN*CP)/( <N/0ENW*CPH) j r6 - H'1ELT/CPW»TMELT # 

R7 ■ HMELT/HTR ; R8 ■ Q»H/KW*TmELT ******************* 

REaO( 5«99S) R1«R2«R3«R4«R5«R6«R8 
R£AD(5«939I NRiMiN 

c *** computational parameters *** 

TIn" 1 «0 
DT-OtOOOS 
ITT-41 
Cm ■ 1.0 
TMELT ■ 1.0 
ITTT ■ 1 
TPhASE ■ R6 
AN ■ N 
AM ■ M 

QLaTM ■ am»an 
AMI ■ AM ♦ 1.0 
Ml ■ M ♦ 1 
M2 ■ Ml + 1 
Nl ■ N + 1 
N2-N+2 
LL»N2¥M1 
Ul ■ 2»M1 + 1 
L • M2 

CD ■ ( R1#¥2 ) » ( AM*»2 ) 

c *** initial values 

SOS ■ 0.0 
SQP ■ 0*0 
SQB ■ 0*0 
SQT"0 . 0 
SQL ■ 0*0 
ULAT • 0*0 
EBOT-q.O 
TTCB-0.0 
MM ■ 0 
TBOT ■ 0*0 

time ■ 0*0 
DO 10 I ■ 1/LL 
F(I) - l.Q 
MELT! n-0 
10 TEXCI I ) - 0.0 
DO 20 J-1/N2 
DO 20 I - 1/Ml 
20 T1 1 1/ J) - TIN 
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TABLE 12. (Continued); 


c ***** horizontal conductances^ vertical COSOUCTANCES 
c /Capacitances ** 
c **** are calculated ******* 
av( 1/1) > o>o 

DO 30 J-2/N1 
OQ 30 1-2/Ml 
QH( 1/ J) - 1.0 
C( 1/ J) - 1.0 
30 GV ( I / J ) - 1.0 
00 40 J ■ 2/Nl 

C(l/J) - (R4»AMl /(a.0*R5*R3» 

GV(1/J) - ( AM»R4 I /( 2.0*R3) 

QH(2/J) - 2*0/(l*0 ♦ AM/(2‘0»R3*R4) ) 

QH(M2/ J)-0.0 
40 GH( 1/ J) - 0.0 
DO 50 I - 2/Ml 
CtI/1) - ( AM'fRl'.^ 4 ) / ( S2*R5 I 
C( 1/N2)-C( I/l ) 

GV(I/2) ■ 2«0/(1.0 + <R1»AM/(R2*R4) ) ) 

QV( I/N2)-GV( 1/2) 

GVl I/l ) - 0.0 
QH( 1/ 1 ) - r4»AM‘/;?1/r2 
50 GH( I/n2)-GH( 1/1 ) 

QHIl/1 ) - 0.0 
GM(l/N2)-0.0 
QH(M2/1 )-0.0 
QH(M2/N2)-0.0 

GHI2/1) ■ 2.0¥R4/( t R2/R1 )*M1«0/(2.0*R3) t l*0/AM)) 
GH(2/N2)-QH(2/1 ) 

GV(1/2) ■ R4/( R3» ( 1 *0/AM ♦ R1/R2)) 

GV(1/N2)-QV(1/2) 

C(l/1) - (Rl¥R4¥( AM¥*2) )/<2*0 ;*r2»R3¥R5) 

C(1/N2)«C(1/1 ) 

DO 55 I-l/Ml 
55 TTC8-TTCB + C( l/l ) 

c **** temperatures at Each node calculated *** 

rtKITE(6/7500) Rl/ R2/ R3/ R4/ R5/ R6/ R8/ 1/ N/ NR 
60 TIME - TIME + DT 
DO 70 I ■ 1/LL 
DO 70 J - 1/Ll 
70 A( 1/ J) -0.0 
DO 80 J-2/N1 
DO 80 I - l/Ml 
K • (J-1)»M1 + I 
A(K/L-M1) - GVl 1/ J) 

AIK/L+Ml) - GVtIi J¥l) 

AlK/L) - -QHlI+l/J) - GVlI/J) - QV(I/J<“l)-<ClI/J)/([)T 
l¥CO) l“GH( 1/ J) 

A( K/L-1 ) - GHl I / J ) 

AlK/L + l ) - QHl 1 + 1/ J) 

80 D(K) - -C( 1/ J ) +T1 1 1/ J)/(CD¥OT) 
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TABLE 12. (Continued) 


DO 90 1-1/Ml 

K-Nl*hl + I „ . . 

A(liL) - -QH(m*l) - - QV1I#2I - C(I/l)/OT 

1*CD» 

A( liL + l ) - QH( I + lil ) 

A(I/L-1) ■ QH( I>1 ) 

A( I/L+Ml ) ■ GV( 1.2) . . 

A(KiL) - -QH(I.N2) - QH(I+1>N2I - GV(1.N2) - CtI*N2l/ 

U0r*CD) 

A( K.L-1 ) ■ GH( I. S2 ) 

A(K.L+1) • QH(I + 1.‘J2) 

90 A » K.L-Ml > ■ QV I I . S2 ) 

DO 110 1-1. Ml 
K-Nl»ni+I 

0(1) - -Cd.ll'fTid.D/lDT'fCOI - R8/(R1»AM) 

110 OK) - -C ( I / N2 ) *Tl I I . N2 ) /( OT*CO » 

0«U ■ -Ct 1. 1 ) *T1 ( 1. 1 I/ICd^OT) - R8/<2*0»^3»=<1 » 

CALL 6ANS0L(A/LLjL1.0*LL) 

DO 120 I • l.LL 
120 T2i I ) ■ D( I ) 

*** QUALITY CALCULATEO *** 


DO l5o J-2.Nl 
DO l5o 1-2. Ml 
K-(J-l)*Ml + I 
innELTK) .EQ*1 ) 30 TO 150 

IF(T2(K) iGT«TMELT)TEXCKI - TEXCIK) + tT2<<» - TMELd 
IF ( T2 K ) « QT • TMELT ) F ( K ) ■ 1*0 “ TExC ( K ) / ( R6*TM£_T ) J T2 
l(K)-TfiELT 

IF(TEXCK).QE*TPhA3E) T2( K )-TMElT + TEXC( < )-T = HASs.<F 
l(K)-O.OihELT(K)-l 


150 continue 

DO 160 J - 1.N2 
DO 160 I - l.Ml 
!< - ( J-1 ) »nl +■ 1 
160 TH 1. J) - T2K) 

QLAT ■ 0.0 
DO 170 J- 2.N1 
DO 170 I - 2.M1 
K - ( J-1 ) '.Ml + I 
170 uLaT - qLAT + 1 1*0 - FK) ) 

IFtAaslULATM - QoAT ) »L r.0.0001 ) GO TO 180 
IF( T1 mE.EQ. 0T ) ITTT ■ ITTT 1 
IF« TIME.EQ.DT ) 30 TO 180 
IF( ITT. NE. ITTT) ITTT - ITTT + 1 
IFl 1 TT.nE.ITTT)G0 to 240 
ITTT - 1 

180 DO 175 I-l.Ml 

175 EBoT-EBOT+1 C( I.D^ITKI.D-HN)) 
TBUT-TIN+(£B0T/TTC3) 

EBUT-O.O 

TRISE - TBOT - TmELT 



TABLE 12. (Continued) 


UCOEF - R8/(TB0T - TMELTl 
WR I TE( 6*8000 ) TIME# TRISE* UCOEF 
c *** energy balanced *** 

DO 200 J - 2#M 

200 SQF - SQF + (Ctl. J)*(Tl(l#JI - TIN)/R6» 

DO 210 I - l#Ml 

SQT-SQT+(C( I#N2I*(T1( I#N2)-TIN)/R6) 

210 SQB - SQB + (Cl I#1 I*(T1' t#l> - TIN)/R6» 

00 220 J ■ 2#N1 
DO 220 I ■ 2#Ml 
K ■ ( j“i ) I 

SQS-SaS+( Cw*( T2( < ) -TIN»*F(K)/R6) 

220 SaL-SQL+(CR‘f(T2(O-TINM‘(l*0-F(K) )/R6) 

OTOT - SQF + SQB + SQS + SQL ♦ QLAT +SOT 

EIn - ( AiH>¥2'^Rl'fR8/’R6HM 1*0 ♦ ( 1 1 0/ ( 2 • 0»R 3 ) ) ) » T I ME 

RATIO - QTOT/EIN 

fmelt-qlat/qlatm 

WRIT£(6#1500)S3F# SOB# SQT# SQS# SQL# QLAT# OTOT# El N# RATIO# 
1 FMELT 
SQF - O.O 
SQB ■ 0*0 
3QT *0 • 0 
SQS ■ 0*0 
SQL ■ 0*0 
TBoT ■ 0*0 

2Aj IF(ABS(QLATM ■ QLAT ) •3T*0 »o 0OI ) GO TO 6o 

998 FORMAT! 7F10«5) 

999 F0RMAT(3I3) 

1000 FORMATdH # 1 3# 1 X # NFl 0 • 5 ) 

1500 FORMATdH #'SQF ■ ' #F10*5#5X# ' SQ8 ■ • # FlO < 5# 5X# • SOT ■' 
1#F10.5#5X# 'SQS ■'#F10*5#5X# 'SOL-dFlO^S#/# '3LAT •' 
2#F10.S# SX# ' QTOT ■ ' # FIO • 5# 5X# ' EIN - ' # FI 0 . S# 5X# ' RAT 10 
3#FiO*5#5X# 'FMElT-' # Fl0»5l 
2000 format (IHI ) 

3000 FORMATdH # I3#NF10.5) 

4000 FORMATdH # T35# ' HORIZONTAL CONDUCTANCE') 

5000 FORMATdH # T35# ' VERTICAL CONDUCTANCE') 

6000 FORMATdH # T3S# ' THERMAL CAPACITANCE') 

7000 FORMAT! IH ) 

7500 • FORMAT ( 1H1#'R1«'#F7«3#3X# ' R2- ' # F6 • 1 # 3X # •R3»'#F3*1#3X 
1# 'R4-'#F7.1#3X# 'R5-'#F5*1#3X# • R6« * # F5 • 3 # 3 X# 'R8"'#F6i3 
2#3X# 'M-'#I2#3X# 'N-'#I3#3X# 'RUN NO# «'#I4) 

8000 FORMAT! IHO# ' TIME • ' # F7 • 5# 5X# ' TR ISE ■ ' # F7 • 5# 5X# • UCOEF 
1 ■' #Fl0i5) 
call EXIT 
END 

Subroutine bansdl!c#n#m#v#io) 
dimension C! ID#M)#V! ID) 
implicit REAL*8! A-H#0-Z) 

L • ( M+1 ) /2 
LI =• L ■ 1 


55 


TABLE 12. (Concluded) 


DO 

2 

IR • 1/Ll 

LH 

■ 

L - IR 

DO 

2 

I ■ 1/LR 

DO 

1 

J ■ 2/ M 

Cl IR/J-1) • Cl IR« J) 


NPl ■ N ♦ 1 - IR 

rtPl ■ M ♦ 1 • 1 

C(NPI#MP1) ■ OtO 

2 Ct IRiM I ■ C(NP1/ *1»l > 

N1 ■ N • I 

00 9 I ■ 1#N1 
IPIV « I 
IRE ■ I + 1 

DO 3 IR > IRE «L 

IF( A0S( Cl IR#1 > ) *wE* ABSlCl IPIV*1 I n 00 TO 3 
IPIV ■ IR 

3 continue 

IF(IPIV.EQ.I) 30 TO 5 

T - vm 

V( I ) • VI IPIV) 

VCIPIV) ■ T 
DO 4 J • 1/M 
T - Cl 1/ J) 

CII/J) - CIIPIV/J) 

4 CIIPIV/J) - T 

5 VII) - VI I )/CI I/l ) 

DO 6 J ■ 2iM 

6 CII/J) - Cl 1/ J)/CI I/l ) 

00 8 IR ■ IRE/L 

T - CIIR/1) 

VI IR) • VI IR) • T»VI I ) 

DO 7 J ■ 2iM 

7 CIIR/J-1) - CI1R*J) ■ T*CCI/J» 

8 Cl IR/M) > 0>0 
IF|L*EQ»N) GO TO 9 
L - L ♦ 1 

9 continue 

VIN) ■ VIN)/CIN/1 ) 

JM - 2 

DO 11 ICE ■ 1/Nl 
IR ■ N - ICE 
DO 10 J ■ 2«JM 
IRfll ■ IR - 1 ♦ J 

10 VIIR) ■ VIIR) • CIIR/J)*VIIRM1» 

IFIJM»EQ>M) GO T3 11 

JM • JH ♦ I 

11 continue 

RETURN 

END 



o r> o o n n 


TABLE 13. COMPUTER PROGRAM FOR BASE EXPOSED TO A STEP 
INCREASE IN TEMPERATURE 


c constant Base teiperature program 

OIHEnSION A|2A3«7)* T1(3/82)« T2(243)« GV(3 

1/82)/ C(3/82)/ F(2A3)/ T£XC(243>/ D<243)/ 'lELT(243) 
IMPLICIT REAL*8 ( a-H/0-2) 

R1 ■ H/H J R2 • H/Sl I R 3 ■ M/S2 ; RA ■ </<»< j ♦»♦**»* 

*** R5 ■ {K/OEN<‘CP)/(<W/OENW*CPW) j r 6 • H**EtT/'CPW»TMELT * 
»#»¥*****♦** 

*»¥ R7 ■ HMELT/HTR I R8 ■ 0*H/KW¥TrtELT ♦¥¥**¥*»♦¥¥¥¥»¥♦¥*¥ 
¥¥¥¥¥¥¥¥¥¥¥¥ 

READ I 5/ 998) R1 / R2/ R3/ R4/ R9/ RS/QELT 
READ(5/999) NR/M/S 

c computational parameters ♦♦♦ 

TIn»i« 0 
OT ■ 0*00050 
ITT-6 
CW " 1*0 
TflELT • 1*0 
ITTT • 1 
TPMASE • R6 
AN ■ N 
AM ■ « 

QLaTM • AM¥AN 
Ml • n ♦ 1 
M2 ■ Ml ♦ 1 
N1 ■ N ¥ 1 

N2*N+2 
LL»N1¥M1 
LI ■ 2»Ml ♦ 1 
L • M 2 

CD • (R1¥¥2)¥< AM»¥2I 

c ¥¥* initial Values 

SQS ■ 0*0 
SQP > 0*0 
SOS • 0*0 

saT**o*o 

SQL ■ 0*0 
QLaT ■ Q*y 
EBOT*o • 0 
TTCB«o *0 

TBOT ■ 0*0 

time ■ 0*0 

OTOTl ■ 0*0 
DO 10 I • 1/LL 
F ( I ) • 1*0 
MELT(I)-0 
10 TEXC(I) ■ OtO 
DO 20 1-1/Ml 
Tl( I/l )-TIN¥DE lT 
DO 20 J-2iN2 
20 Til 1/ J) « TIN 


57 



TABLE 13. (Continued) 


C HORIZONTAL CO'^OJCTA.xCES* VERTICAL CONDUCTANCES 

C ^CAPACITANCES ** 

C are CALCULATED **♦*»*♦ 

QVdil ) ■ 0«0 
00 3 O J"2*N1 
OO 30 I«2#hl 
QH( li J) ■ 1.0 
C ( 1 4 J ) ■ 1.0 
3 O QV ( I i J ) ■ 1.0 
00 40 J - 2iNl 

C(l/J» • (R4¥AHl/(2.0»R5^fR3l 
QVtl^J) ■ (AM*R4)/(2*0'»R3) 

GM(2/J) ■ 2*0/(1.0 + Ah/( 2*0»R3¥R4» I I 
QM(M2> j)«OiO 
4 O QH(l.J) ■ 0*0 
00 50 1 - 2,fil 
C(I#l) ■ ( Ah»Rl»R4 ) /( R2*R5) 

Cl I#N2)»C( 1#1 I 

□ Vlli2) ■ 2»0/(1.0 + (R1*A(1/(R2*R4) >) 

QVl I/N2)«GVl 1*2) 

QV( 1*1 ) ■ 0.0 
GH(1*1) ■ R4»AH*R1/R2 
50 GHl I*N2)-QH( 1*1 > 

GH (1*1) ■ 0 . 0 
QH11*n2)-0.0 
GH(H2*1 )«0*0 
QM(H2*N2)«0*0 

GH(2*1) ■ 2*0»R4/l t R2/R1 )*( 1'0/(2*0*R3) ♦ 1*0/AM)) 
QHl2*N2)«iQH(2*l ) 

GV(1*2) - R4/(R3j*(1.0/AH ♦ R1/R2)) 

QV(1*N2)«GV(1*2) 

Cll*l) ■ (R1»R4»( AM*»2) )/(2»0*R2*R3*R5) 

Cll*N2)-C( 1*1 ) 

OO 55 I«l*ril 
55 TTCB-TTCB + Cl 1*1 ) 

C **** temperatures at each node CALCULATED 
WRITE (6*7250) NR 

wRITE(6*7500) R1* R2* R3* R4* R5* R6* DELT* 1* N 
60 TIME - TIME + dT 
00 70 I - l*LL 
00 70 J - 1*L1 
70 A( I* J) -0.0 
00 80 J-2*Nl 
00 80 I ■ 1,M1 
K - (J-2)»M1 ♦ I 
IF (J.NE»2) GO TO 75 

OlK) - -C( I* J) i>Tl ( I* J)/( Cd»0T)-3V( I* J ) *T1 ( I* J-1 ) 

GO TO 79 

75 DIK)- -C(I*J) ♦ T1 ( I* J)/(CO*OT| 

A(K*L-P11 ) - QV( 1* J) 

79 A(K*L + ni ) ■ QV( r* J-t -1 ) 
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TABLE 13. (Continued) 


AlrCiL) - -uHlI + liJI - GV(I,J) - OVlI»Jn)-lCtI,J)/(DT 

A U< L-1 ) ■ GH( 1/ J ) 

ao A( K/ L + a ) • QH t I + l i J ) 

UO 90 
K« N»fll+-I 

AlK;L) - -QH(I*S21 - GH<1^1 >n2) - GVII/nS) - ClIiN2)/ 

1 ( OT»CO ) 

Al K jL-1 ) ■ GH( I * n2 ) 

AIK/L + 1 ) « GH( I + 1»n21 
30 A 1 K> L-ril ) - GV( 1 » N2 » 

DO 110 I«l>f11 
<■ r-j^ni + I 

110 DU) - -Cl UN2)¥T1 ( I/N2t/<uT»CD» 

CALL BANSUL ( A * LO i L 1 * D * LL ) 

DO 120 I - 1#LL 
120 rati) • Dill 
*** duality calculated *** 

00 150 J»2»Nl 


DO 150 1-2/Ml 

<■ ( J-2 ) ♦Ml + I 

IFImElTU) /EQ*! I 30 TO ISO 

IF(T2tK)«GT*TnELT)TEKCUl - 
IF(T2(K) .QT/TriELT)FU) - 1‘0 


TExCl K ) ♦ ( T2( <1 - TMELT I 
- TEXCU)/1=16*TME_TI < T2 


^IFl TEXCU) •GE.TPMA3EIT2(Kl“TnELT+T£XCU)-T^HASE)F 

hk)-o.o;melt(ki»i 
1 fj 0 C 0 :■< T l N U E 

DO 160 J - 2/N2 
DO 160 I - 1/Ml 
< - ( J-2 ) »Ml I 
16'J TKI/JI - T2U) 

QLaT » 0*0 
DO 170 J- 2/Nl 
DU 170 I - 2/Ml 
K ■ (J-2)^M1 + I 
I7u OLaT « QLaT + ( 1*0 “ FlOl 

T8QT - 0*0 
DO 175 I-l/Ml 

175 EBjT“£BOT+lC(l/lU<Tl(l<ll*I’lNl • 

TBJT-TIN+ I EB0T/TTC3 ) 

E8uT-0 *0 

TRISE - T80T - TMELT 
♦♦♦ energy balanced ♦♦♦ 

SOF - 0.0 
SQB - 0*0 
SQT ■ 0*0 
SQS ■ 0*0 
SQL - 0‘0 
DO 2u0 J » 2/ Nl 

200 SOF ■ SQF + lCtl/JUtTl(l/J> - TIN)/R6) 



TABLE 13. (Continued) 


00 210 I ■ liMl 

SQT»SaT+(C( I#N2)¥(T1( 1 * N2 I -TIN ) /R6 I 
210 SQB • 8QB ♦ ( C ( 1 • 1 > » < T1 1 1 < H • TIN)/R6> 

00 220 J ■ 2iNl 
00 220 I ■ 2iMl 
K • (J-2)¥M1 + I 

SQS-SQS+ (CW¥(T2(0-TIN)<>F(K)/R6) 

220 SQL-SQL+ ( CW»( T2( <) -TI N » • ( 1 »0-F ( K 1 )/R6) 

QToT2 ■ SQF + SQ3 ♦ SQ8 ♦ SQL ♦ QLAT ♦SQT 


QRATE - (QT0T2-QT0T1 )/DT 

UCOEF • (QRATE#RS)/( (TB0T-TMELT»*(M*»2J*Rl» 
l/(2*R3) I 
QTOTI-QT0T2 
FMELT ■ QLAT/(AM*AN) 

IF (ABSIQLATM - 3wAT ) •LT'O'OOOl 1 QO TO 180 
1TIME»EQ»DT) ITTT ■ ITTT ♦ 1 
1TIME*EQ*DT) 30 TO 180 
( ITT»NE»ITTT) ITTT • ITTT 
30 TO 2 a0 


(1+2*R3» 


♦ 1 


IF 
IF 
IF 

IF ( ITT.NE»ITTT) 

ITTT - 1 

write (6#8000) TIlE/TRlSEiUCOEF t 

WRITE <6<1500) SaF»8aB/8QT>SaS/SQL#QLAT*aT0T2#FMELT 
IFtABStQLATH - QLAT ) • 3T»0t0001 > GO TO 60 
FORrlATIlH /lOXi'SQF ■ ' # FlO • 5X* * 8Q8 -’^FlO^B/BX* SOT 

1 -SFlOiS/BX# /<11X/'SQS -'^FlO^S/SXi 'SOL -’^FlO'BMX 
2#'QLAT -SFlOiSz/i lOx# *QTOT •• /F10»S/ 4X# ‘ FIELT • 
3#F10*5I 

998 F0RMATI7F10*5) 

F0rMAT( 3I31 
FORMATdHl I 

FORMAT (IHlMXi’R'JN NO* S I4) 

F0RMAT(1H0/9X# *R1 ■ * i P7 • 3# 3Xi ' R2 ••#F6*1*3X#R3» 
l»F5*l* 3X# ' R4 •' /F7*l#/*10Xi 'R5 -» #F5»1< 5X« ' R6 »dF5»3 
2#4X/*0ELT -dFe'^^aXi'M - '• I2/5X*'N -'<I4» 

BOOO FORMAT! 1H0#9X/ ' time ■ dF7»5*5Xi ' TRIBE ■•#F7»5#3X 
1#'UC0EF -'iFlC.BI 
WRITE l6/2000> 
call EXIT 
END 

SUBROUTINE BANSOL ( C* NNEW* MNEw< V< ID* 

IMPLICIT REAL»8( A-H/O-Z J 
dimension C(ID/MI*V(ID) 

L - (MNEW+1 1/2 
Ll ■ L - 1 


180 

240 

1500 


999 

2000 

7250 

7500 


DO a IR ■ l*Ll 
LR • L - IR 
DO 2 I ■ 1#LR 
00 1 J • 2^f1NEW 
1 C(IRO-l) - C( IR» J) 
NPl ■ NNCW + 1 - IR 
MPl • MNEW +1-1 
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table 13. (Concluded) 


C(NP1>MP1) « 0»Q 
a C( IR/MNEW) • C(NPl»rtPl) 

N1 - NNEW ■ 1 
00 9 I - 1*N1 
IPIV ■ I 
IRE - I + 1 
00 3 IR - IRE <L 

IF( ABS(C( IRi I ) ) •t.E'ABSlCi IPlV«l ) ) > QO TO 3 
IPIV - IR 

3 continue 

IF( IPIVtEQ.l ) GO TO 5 
T • V( I ) 

V( I ) ■ V( IPIVI 
V( IPIV I • T 
00 4 J - liHNEW 
T - C(I/J) 

C(I*J) - CdPiV* J1 

4 C( IPIVi J) ■ T 

5 V<I) - Vdl/Cdd) 

00 6 J ■ SiNNEW 

6 Cd/ Jl ■ cd# J)/Cd#l > 

DO 8 IR ■ IRE/L 

T - CdRill 

VdRI ■ VdR) • T»V(I) 

DO 7 J ■ 2>NNEW 

7 CdR/J-ll - CdR/JJ - T*Cd*J) 

8 CdR^MNEW) • 0*0 

iFt L*EQ«NNEW ) GO TO 9 
L ■ L 1 

9 CONTINUE 

V(NNEW) ■ VINNEW)/C(NNEWiH 
Jpl • 2 

DO 11 ICE • liNl 
IR - NNEH ■ ICE 
00 10 J ■ ztjn 
IRMl • IR • 1 + J 

10 VdR) ■ VdRI - CdR* J)*VdRMl ) 

IF( JrttEQ*MNEW) QO TO 11 

JM ■ JM + 1 

11 continue 

RETURN 

END 
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IF(T2(K).LT.TMELT)TEXC(K)=TEXC(K)-(T2(K)-TMELT) 

1F(T2(K).LT.TMELT)F(K)=1.0-TEXC(K)/(R6*TMELT);T2(K)=TMELT 

IF(TEXC(K).GE.TPHASE)T2(K)=TMELT-TEXC(K)+TPHASE;F(K)=0.0;MELT(K) 

= 1.0. 

c. In statement 170, change sign to give QLAT = QLAT - (1.0 - F(K)). 

d. In the line above statement 220 and statement 220, make changes to give 
SQS=SQS+(CW*(T2(K)-T1N)*(1.0-F(K)/R6) 

220 SQL=SQL+(CW*(T2(K)-TIN)*F(K)/R6) 


An extensive parametric study for the case of melting with the initial PCM 
temperature at the fusion temperature was conducted using the computer program given in 
Table 12. This parametric study is discussed in more detail and the results are presented in 
Section IV. 


IV. GENERAL PARAMETRIC STUDY 


A. Introduction 

Presented in this section arc some details and typical results of a parametric study 
for a particular phase change device application. Since there are many phase change 
device applications, each having its own unique constraints, it would be unrealistic to 
propose that the results obtained for a particular set of conditions and constraints could 
be quantitatively extrapolated to other conditions and constraints. However, results 
obtained for a particular case such as those presented herein may afford valuable insight, 
qualitatively, into such considerations as trends and acceptable operating ranges. 

) The parametric study consisted of a thermal analysis of the symmetrical half-cell 
shown in Figure 29 with the base exposed to a constant heat flux density. The model 
and assumptions are described in Section III. 


B. Scope and Definitions 

Numerical results liave been obtained using the implicit formulation described in 
Section III. 
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A band algorithm [14] was used to reduce the required computer storage m 
comparison to that which would have been required using a general matrix inversion 

routine. 

Of the seven dimensionless variables previously given, some of the physical 
parameters were considered to be fixed for the computations. Independent variation of 
all seven parameters would render the problem untractable. 

By selection of a particular PCM and metal for the housing and fins, the three 
dimensionless parameters denoted by R 4 , R 5 . and R^ are all fixed. Since aluminum was 
typically used in most of the previous phase change device applications and paraftin was 
commonly used as the PCM, properties corresponding to these materials were used in t e 
computations. For the results presented herein, the representative properties used in the 
computations for these materials are given in Table 14. With the property values given in 
Table 14, the values for the fixed dimensionless ratios are 


R 4 = 1069 

Rs = 753 


and 


R. = 0.381 


Values of the geometrical variables and the heat flux density chosen for the numerical 
predictions are based on a survey of past applications. The values are listed in Table 15. 
For the range of variables outlined in Table 15, there are 300 possible c-ombinations 
which indicates the need for an equal number of computer runs to obtain results 
encompassing these conditions. Because of the large number of computer runs involved it 
was not practical to conduct a sensitivity study for each run on the choice of M and N 
and DT, a suggestion outlined in Section III under the instructions for using he 
computer programs. The choice of M and N was made such that the node size S for the 
range of geometric variables given in Table 15 never exceeded an approximate value of 
0.254 cm (0.1 in.). This value, for some cases, was not sufficwntly small to avoid there 
being an effect on the character of the transient results. This is discussed later m this 

section. 

Specifically, for the results shown hereafter in this section, the parametric 
computations were’ performed for the case of melting only. The entire cell was considered 
to be initially at the fusion temperature of the PCM and in the solid phase. 

If one assumes that the phase change is reversible, the results can also be applied 
to freezing There is evidence, however, that the rate ot freezing is typical y 
underpredicted by the model. However, in most designs, tliese data will indicate 
conservatively slower freeze rates than actually occur. Additional discussion of the 
freezing case is presented in Section V. 
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TABLE 14. PCM AND METAL PROPERTIES USED IN COMPUTATIONS 



Den 

sity 

Spe 

He 

cific 

at 

Thermal 

Conductivity 

Heat of 
Fusion 

Melt 

TemDerature 

Material 


/lbm\ 

\w) 

kJ 

/ Btu \ 

J 

/ Btu \ 

kJ 

7Btu\ 


(R) 



kg-°C 


sec-m-°C 

\hr-ft-°F/ 

kg 

\ Ibm/ 

K 

Metal 

2707 

(169) 

0.88 

(0.21) 

160.9 

(93) 

- 


— 

_ 

PCM 

80 

( 50) 

2.09 

(0.50) 

0.15 

(0.087) 

232.4 

(100) 

291.7 

(525) 


TABLE 15. VALUES OF GEOMETRIC VARIABLES AND HEAT FLUX DENSITY 
CHOSEN FOR THE NUMERICAL STUDY 


Cell Height, H 
One Half-Fin Spacing, W 
Fin Thickness, S 2 
Base Plate Thickness, S, 
Heat Flux Density, q" 


cm 

1.27 

2.54 

(in.) 

(0.5) 

(1.0) 

cm 

0.127 

0.254 

(in.) 

(0.05) 

(0.1) 

cm 

0.0127 

0.0508 

(in.) 

(0.005) 

(0.020) 

cm 

0.0635 


(in.) 

(0.025) 


J 

s-m^ 

1076.4 

2152.8 

Btu 

h-ft^ 

(341.4) 

(682.9) 


5.08 

10.16 

20.32 

(2.0) 

(4.0) 

(8.0) 

0.508 

1.016 

2.032 

(0.2) 

(0.4) 

(0.8) 

0.127 



(0.050) 




4305.6 6458.3 

(1365.7) (2048.6) 


The problem is transient in character. The computer program provides a printout 
of a number of computed quantities at time intervals which can be chosen as any 
multiple of the time increment used in the numerical expressions. The quantities 
computed at each time included the temperature difference between the base and the 
melt temperature (ATj,), an overall heat transfer coefficient (U), the latent energy which 

represents the amount of PCM which has melted (QLAT), and the sensible energy stored 
in the fin (SQS), the base (SQB), the cover (SQT), the solid PCM (SQS), and the liquid 
PCM (SQL). Computations automatical'y terminate when all the PCM has melted. 

For the geometry considered here, the latent energy associated with the PCM per 
unit area of the base is expressed by 
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( 22 ) 



2 p hf H 
Y+ S2/W 


This parameter is defined because it is useful in expressing the latent energy storage 
capability of the device. In using equation (22) to determine the latent energy storage 
capability, the appropriate area to use with this equation is the inside plan area of the 
device. This leads to a negligible error since equation (22) is based on a symmetrical cell 
and the two outermost end cells of the device are unsymmetrical. The latent energy is 
linearly related to H for fixed values of Sj and W. 


C. Results 

For the range of H, W, S 2 , and the physical properties considered in this study, 
equation (22) is represented graphically in Figures 33 through 35. Also shown is the case 
which corresponds to the case for which the entire housing is filled with 
PCM without fins. The effect of increasing fin thickness is readily apparent by comparing 
the curves shown in Figures 33 through 35. 

An important consideration, especially in temperature control design, is the 
temperature rise of the base. Ratios of the base temperature rise at the termination of 
melting to the absolute fusion temperature of the PCM for the range of conditions 
considered herein are shown plotted versus W in Figures 36 through 50. Over the 
examined range, the base temperature rise increases with an increase in W. For a specified 
set of parameters, the data indicate that the temperature rise approaches an asymptotic 
value which corresponds to expected behavior since the problem should reduce to that of 
a one-dimensional case as W becomes large. It is also expected that the temperature rise 
would increase at very small values of W as the problem approaches that for the heating 
of a metal slab. Computations were not performed for these very low values of W. It 
should also be noted that some of the values for the base temperature rise which are 
shown are too large to be of much practical significance; some of these are included, 
however, to indicate trends and to facilitate approximate interpolation. For some of the 
larger heat-flux-density cases, the curves are terminated at lower values of W than they 
are for the lower heat-flux-density cases because the resultant temperature rises are much 
too large for practical considerations. 

In Figures 36 through 50 as well as in all subsequent figures showing results of 
the parametric study, computed data are represented by solid circles. Curves have been 
faired through the points; it should be recognized, however, that the value for points on 
the curve between two computed values are only approximate. 


The data shown in Figures 36 through 50 are cross plotted versus fin height in 
Figures 51 through 62. For some geometries within the considered range, the curves 
indicate the presence of a relative maximum and a relative minimum. This feature is 
aUributable to the two-dimensional character of the problem. This characteristic is 
discussed in more detail in Appendix B. 
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Figure 36. Base temperature rise at the termination of melting versus W for 
H = 1 .27 cm (0.5 in.) and 83 =0.1 27 cm (0.05 in.). 
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Figure 37. Base temperature rise of the termination at melting versus W for 
H = 2.54 cm ( 1 .0 in.) and S 2 =0.1 27 cm (0.05 in.). 



5.08 cm (2.0 in.) 
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Figure 38. Base temperature rise at the termination of melting versus W for 
H = 5.08 cm (2.0 in.) and Sj = 0.127 cm (0.05 in.). 
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Figure 39. Base temperature rise at the termination of melting versus W for 
H = 10.16 cm (4.0 in.) and S 2 =0.127 cm (0.05 in.). 
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Figure 40. Base temperature rise at the termination of melting versus W for 
H = 20.32 cm (8.0 in.) and Sj =0.127 cm (0.05 in.). 




> 1.27 cm (05 in.) 

- 0.0508 cm (0.02 in.) 6458.3 (2048.6) 
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Figure 41 . Base temperature rise at the termination of melting versus W for 
H = 1.27 cm (0.5 in.) and S 2 = 0.0508 cm (0.02 in.). 




75 


Figure 42. Base temperature rise at temiination of melting versus W for 
H = 2.54 cm (1.0 in.) and S 2 = 0.0508 cm (0.02 in.). 
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Figure 44. Base temperature rise at termination of melting versus W for 
H = 10. 16 cm (4.0 in.) and Sj = 0.0508 cm (0.02 in.). 


6458.3 (2048.6) 



;ure 45. case temperature rise at termination oi melting versu! 
H = 20.32 cm (8.0 in.) and S 2 = 0.0508 cm (0.02 in.). 



'1.27 cm (0.5 in.) 

~ 0.0127 cm (0.005 in.) 6458.3 
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Figure 46. Base temperature rise at termination of melting versus W for 
H= 1.27 cm (0.5 in.) and = 0.0127 cm (0.005 in.). 
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Figure 47. Base temperature rise at termination of melting versus W for 
H = 2.54 cm (1.0 in.) and S 2 = 0.0127 cm (0.005 in.). 
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Figure 48. Base temperature rise at termination of melting versus W for 
H = 5.08 cm (2.0 in.) and = 0.0127 cm (0.005 in.). 




10.16 cm (4.0 in.) 
0.0127 cm (0.005 in.) 
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Figure 49. Base temperature rise at termination ot melting versus W for 
H= 10.16 cm (4.0 in.) and Sj =0.0127 cm (0.005 in.). 



20.32 



Figure 50. Base temperature rise at termination of melting versus W for 
H = 20.32 cm (8.0 in.) and S 2 = 0.0127 cm (0.005 in.). 






1076.4 J/im^ (341.4 Btu/h-ft^) 
0.0127 cm (0.005 in.) 
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Figure 5 1 . Base temperature rise at the termination of melting versus H for q" = 1076.4 
J/s-m2 (341.4 Btu/h-ft^ ) and Sj = 0.0127 cm (0.005 in.). 



1076.4 Hfttr (341.4 Bta/h-ft^) 



Figure S2. Base temperature rise at the termination of melting versus H for q 
J/s-m* (341.4 Btu/h-ft*) and Sj = 0.0S08 cm (0.020 in.). 





= 1076.4 (341.4 Btu/h-ft 

= 0.127 cm (0.050 in.) 



gure 53. Base temperature rise at tlie termination of melting versus H for q" = 1 076.4 
J/'s-m^ (341.4 Btu/h-ft^ ) and Sj =0.127 cm (6.050 in.). 
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Figure 54. Base temperature rise at the termination of melting versus H for q" = 2 152.8 
J/s-m^ (682.9 Btu/h-ft^) and Sj = 0.0127 cm (0.005 in.). 





- 21 52.8 (682.9 Btu/h-fT ) W cm (in.) 

« 0.0508 cm (0.020 in.) ■ 
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Figure 55. Base temperature rise at the termination of nielting versus H for q" = 2152.8 
J/s-m* (682.9 Btu/h-ft* ) and Sj = 0.0508 cm (0.020 in.). 



215Z8 (682J Btu/h^) 



Figure 56. Base temperature rise at the termination of melting versus H for q 
J/s-m^ (682.9 Btu/h-ft^ ) and Sj =0.127 cm (6.05 in.). 



4305.6 (1365.7 Btu/h-ft^) W cm (in.) 

0.0127 cm (0.005 in.) I 
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Hgure 57. Base temperature rise at the termination of melting versus H for q 
J/s-m* (1365.7 Btu/h-ft^ ) and Sj = 0.0127 cm (0.005 in.). 



= 4305.6 Mi-vcr (1365.7 Btu/h ft'^) W cm (in.) 

= 0.0508 cm (0.020 in.) I 
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Figure 58. Base temperature rise at the termination of melting versus H for q" = 4305.6 
J/s-m^ (1365.7 Btu/h-ft^) and S 2 = 0.0508 cm (0.020 in.). 





4305.6 (1365.7 Btu/h-ft^) 






6458.3 ilt-nr (2048.6 Btu/h-ft^) 
0.0127 cm (0.005 in.) 
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Figure 60. Base temperature rise at the termination of melting versus H for q 
J/s-m^ (2048.6 Btu/h-ft^ ) and Sj = 0.0127 cm (0.005 in.). 



= 6458.3 (2048.6 Btu/h ft 

- 0.0508 cm (0.020 in.) 
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Figure 6 1 . Base temperature rise at the termination of melting versus H for q" = 6458.3 
J/s-m^ (2048.6 Btu/h-ft^ ) and Sj =0.0508 cm (O.O20in.). 




6458.3 (2048.6 

0.127 cm (0.050 in.) 
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Figure 62. Base temperature rise at the termination of melting versus H for q “ 6458.3 
J/s-m^ (2048.6 Btu/h-ft^) and Sj =0.127 cm (0.050 in.). 



The base temperature rises shown in Figures 36 through 62 correspond to the 
termination of melting. If the entire amount of PCM is not melted, the temperature rise 
of the base may be considerably less than that reflected in Figures 36 through 62. As an 
indication of the transient behavior of the base temperature, certain transient data are 
plotted in Figures 63 through 77. For some cases, the base temperature rise remains 
relatively low over a large percentage of the phase change process and rises rapidly near 
the end of the process. In other cases, the rise occurs nearly uniformly throughout the 
phase change process. In either case, designs based on the data shown in Figures 36 
through 62 should be conservative, especially if the design does not require complete 
melting of the PCM. However, as indicated by these data, in cases where base plate 
temperature control is the object, it may be more desireable to limit phase change 
in order to avert unreasonable fin spacing and/or heights. 


Some of the curves shown in Figures 63 through 77 exhibit fluctuations in the 
transient base temperature. For those curves where a sharp rise occurs in the computed 
base temperature (e.g., the upper curve in Figure 63), the accuracy with which the 
predictions describe the actual physical variation is questionable, particularly during the 
ear y part of the transient. The computing model attempts to predict the behavior of a 
continuous system with use of a limited number of finite size nodes. Since the phase 
change is handled in the model by forcing an entire node to remain at the phase change 
temperature until sufficient energy has accumulated to account for the phase change one 
would expect some stepping in the results unless the nodes are made very small. This 
however, is related to the choice of M and N, which in turn dictates computed 
requirements. Since the problem examined in the parametric study involved a step input 
of a constant heat flux at the base, the initial level to which the computer temperature 
rises IS dependent on the node size. The influence of N on the computed results is 
illustrated in Figure 63 for three runs, in Figure 64 for one run, and in Figure 65 for one 
run. For these selected runs, computations are shown for N values of 5 and 10. When 5 
IS used, a horizontal row of nodes represents 20 percent of the PCM. Obviously the 
precise detail of the transient variation is affected in each case, but the principal 
difference is at the beginning; it appears that increasing N results in a smoothing out of 
the fluctuations and shortens the initial plateau. The general trend of the results remains 
consistent and there appears to be only a slight effect in the latter stages of a run. For 

reference purposes the values of N used in the computations are given for the respective 
curve in Figures 63 through 77. ^ 


.n graphically in Figures 78 

through 89. These times exceed those which would be required by purely latent 
considerations because of the effect of sensible heating. 

Since there is a certain amount of sensible energy storage, ratios of the total 
stored energy at the termination of melting to the latent energy required to change phase 
shown plotted versus heat flux density at the base in Figures 90 through 


A few additional computer runs were made for selected values of the parameters 
to examine the influence of the fusion temperature on the predicted temperature rise of 
the base at the termination of melting. The results of these computations are shown in 
Figures 105 through 107. 
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Figure 63. Transient base temperature rise versus fraction of PCM melted 
for H = 1.27 cm (0.5 in.) and Sj = 0.0127 cm (0.005 in.). 





Figure 64. Transient base temperature rise versus fraction of PCM melted 
for H = 1.27 cm (0.5 in.) and Sj = 0.0508 cm (0.020 in.). 
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Figure 66. Transient base temperature rise versus fraction of PCM melted 
for H = 2.54 cm (1 in.) and Sj = 0.0127 cm (0.005 in.). 
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Figure 68. Transient base temperature rise versus fraction of PCM melted 
for H = 2.54 cm (1 in.) and S 2 =0.127 cm (0.05 in.). 










.6 0 


FRACTION MELTED 


Figure 70. Transient base temperature rise versus fraction of PCM melted 
for H = 5.08 cm (2 in.) and Sj = 0.0508 cm (0.02 in.). 
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Figure 7 1 . Transient base temperature rise versus fraction of PCM melted 
for H = 5.08 cm (2 in.) and Sj =0.127 cm (0.05 in.). 
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Figure 73. Transient base temperature rise versus fraction of PCM melted 
for H = 10.16 cm (4 in.) and = 0.0508 cm (0.02 in.). 
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Figure 75. Transient base temperature rise versus fraction of PCM melted 
for H = 20.32 cm (8 in.) and Sj = 0.0127 cm (0.005 in.). 
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Figure 77. Transient base temperature rise versus fraction of PCM melted 
for H = 20.32 cm (8 in.) and Sj =0.127 cm (0.05 in.). 
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q" = 1076.4 J/i-m^ (341.4 Btu/h-ft^) 
$2 ^ 0.0508 cm (0.02 in.) 



Figure 79. Time required for PCM to melt for q” - 1076.4 J/s-m^ 
(341.4 Btu/h-ft^ ) and Sj = 0.0508 cm (0.02 in.). 
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q" = 21 52.8 J/s-m^ (682.9 Btu/h-ft^) / 

So = 0.0508 cm (0.02 in.) / -S' 

^ /' 5 * 



01 23456789 


H, in. 

Figure 82. Time required for PCM to melt for q"' = 2152.8 J/s-m^ 
(682.9 Btu/h-ft^ ) and = 0.0508 cm (0.02 in.). 
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Figure 83. Time required for PCM to melt for q" - 2152.8 J/s-m^ 
(682.9 Btu^-ft^ ) and =0.127 cm (0.05 in.). 
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Figure 84. Time required for PCM to melt for q” = 4305.6 J/s-m 
(1365.7 Btu/h-ft^) and Sj = 0.0127 cm (0.005 in.). 





















H, in. 

Figure 87. Time required for PCM to melt for q ' = 6458.3 J/s-m" 
(2048.6 Btu/h-ft^ ) and Sj = 0.0127 cm (0.005 in.). 
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Figure 90. Fnergy ratio versus heat flux density for W = 0. 127 cm 
(0.05 in.) and Sj = 0.0127 cm (0.005 in.). 
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Figure 9 1 . Energy ratio versus heat flux density for W = 0. 1 27 cm 
(0.05 in.) and Sj = 0.0508 cm (0.02 in.). 
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Figure 93. Energy ratio versus heat flux density for W - 0.254 cm 
(0.1 in.) and Sj = 0.0127 cm (0.005 in.). 
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Figure 95. Fnergy ratio versus heat tlux density tor W = 0.254 cjn 
(0. 1 in.) and S, =0.1 27 ein (0.05 in.). 
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Figure 96. Energy ratio versus heat flux density for W = 0.508 cm 
(0.2 in.) and S 2 = 0.0127 cm (0.005 in.). 
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Figure 97. Energy ratio versus heat llux density for W = 0.508 cm 
(0.2 in.) and S 2 = 0.0508 cm (0.02 in.). 
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Figure 99. Energy ratio versus heat flux density for W = 1.016 cm 
(0.4 in.) and = 0.0127 cm (0.005 in.). 
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Figure 100. Energy ratio versus heat flux density for W = 1.016 cm 
(0.4 in.) and Sj = 0.0508 cm (0.02 in.). 
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Figure 101. Energy ratio versus heat flux density for W = 1 .0 1 6 cm 
(0.4 in.) and S 2 =0.127 cm (0.05 in.). 
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Figure 103. Energy ratio versus heat flux density for W = 2.032 cm 
(0.8 in.) and Sj = 0.0508 cm (0.02 in.). 


137 








^ «o ‘liv 


E 

U 

5 


-l“ 
o 


^/5 


(O 


o 


LT) 

d 



CO 

d 


CVI 

d 


139 


Figure 105. Base temperature rise versus W for three different fusion temperatu 
and H = 1 .27 cm (0.5 in.) and Sj = 0. 1 27 cm (0.05 in.). 
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Figure 106. Base temperature rise versus W for three different fusion temperatures 
and H = 2.54 cm ( 1 .0 in.) and S 2 =0.1 27 cm (0.05 in.). 



0.20 
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Figure 107. Base temperature rise versus W for three different fusion temperatures 
and H = 5.08 cm (2.0 in.) and Si =0.127 cm (0.05 in.). 



Some representative examples illustrating how the results can be used by a 
designer arc given in the following section. For conditions significantly different from 
those used in the computations (Tables 14 and 15), it is suggested that the program 
described in Section 111 be used to determine the required design information. 


D. Examples 

Some hypothetical examples which a thermal designer may encounter illustrating 
use of the results are outlined as follows. All examples assume the designer has selected a 
housing design and materials as discussed earlier. 

Exajnple 1 

Given: The heat flux density is specified to be 1076.4 J/s-m^ (341.4 Btu/h-ft^) 
and the designer wishes to limit the temperature increase to 5.56 K (10°R). The PCM 
melt temperature is 291.67 K (525°R), using these inputs the maximum permissable base 
temperature rise ratio is computed to be 


_ 5.56 

~ 291.67 


19.05 X 10‘3 


Case a The designer wishes to use a fin thickness of 0.0127 cm (0.005 in.). 

• Find an acceptable geometry to satisfy the given conditions. From Figure 
51, a height of 2.54 cm (1 in.) can be used with a half-fin spacing of 
0.127 cm (0.05 in.) to achieve the desired conditions. 

• Determine the energy stored. From equation (22) or from Figure 33, the 
latent energy storage is 4506 kJ/m^ (397 Btu/ft^). The results shown in 
Figure 90 show that the total energy stored at the termination of melting 
for these conditions is approximately 5 percent higher than the latent 
energy stored, the difference being sensible energy storage. 

Case b - The designer wishes to use a fin thickness of 0.127 cm (0.05 in.). 

• Find an acceptable geometry. From Figure 53, it is seen that a fin height 
of 10.16 cm (4 in.) with a half-fin spacing of 0.254 cm (0.1 in.) can be 
used. 

• Determine the energy stored. From equation (22) or from Figure 35, the 
corresponding latent energy storage is 15 128 kJ/m^ (1333 Btu/ft^). At 
the termination of melting, the toial stored energy exceeds the latent 
value by approximately 5 percent as shown in Figure 95. 
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Example 2 

Given: The designer is limited by space to a plan area for the phase change device 
of 0.0697 m^ (0.75 ft^); the heat flux density is specified to be 2152.8 J/s-m^ (682.9 
Btu/h“ft^); the duty cycle of the controlled medium requires latent energy 948.9 kJ (900 
Btu) of storage. For these conditions, latent energy storage per unit area is 

948.9 

E^: = = 1.3614 X 10^ kJ/m^ (1200 Btu/ft^) 

K 0.0697 

Case a — The designer wishes to use a fin thickness of 0.0127 cm (0.005 in.). 

• Find the geometry. From Figure 33, it is noted that the required height is 
approximately 7.62 cm (3 in.) regardless of fin spacing over the range 
shown. 

• Determine the base temperature rise. Assume a half-fin spacing of 0.254 
cm (0.1 in.) and a height of 7.62 cm (3 in.). From Figure 54, the 
corresponding base temperature rise at the termination of melting is 


AT^ ^ 0.132 X 291.67 ^ 38.5K(69°R) 

Case b — The designer selects a fin thickness of 0.127 cm (0.05 in.). 

• Find the geometry. From Figure 35, it can be seen that there is a wider 
range of fin height-spacing combinations which will satisfy the latent 
energy storage requirements. 

• Determine the base temperature rise. First, assume a fin height of 10.92 
cm (4.3 in.) and a half-fin spacing of 0.127 cm (0.05 in.). From Figure 
56, the corresponding temperature rise of the base at the termination of 
melting is approximately 

AT^ ^ 0.027 X 291.67 ^ 7.88K(14°R) 


• Determine the base temperature rise for a second acceptable geometry. 
Second, assume a fin height of 8.38 cm (3.3 in.) and a half-fin spacing of 
0.508 cm (0.2 in.). From Figure 56, the corresponding temperature rise of 
the base at the termination of melting is approximately 


AT^^ ^ 0.06 X 291.67 ^ 17.5K(32°R) 
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Example 3 


Given: In a phase change device the designer wishes to utilize 0.0127 cm (0.005 
in.) thick fins and he wants to limit the base temperature rise to 13.89 K (25‘'R). 
Consequently, the base temperature rise ratio is not to exceed 



13.89 

291.67 


= 48 X 10'3 


Case a A heat flux density of 1076.4 J/s-m^ (341.4 Btu/h-ft^) is required. 

• Determine acceptable geometry and the corresponding latent energy 

storage capacity. From Figure 46, a height of 1.27 cm (0.5 in.) can be 
used with a half-fin spacing up to approximately 0.48 cm (0.19 in.). The 
corresponding latent energy storage capacity, from equation (22), is 2338 
kJ/m^ (206 Btu/ft^). 

• Determine acceptable geometry and the corresponding latent energy 

storage capacity. From Figure 51, heights up to 3.81 cm (1.5 in.) for a 
half-fin spacing of 0.254 cm (0.1 in.) and up to 9.65 cm (3.8 in.) can be 
used for a half-fin spacing 0.127 cm (0.05 in.). The corresponding upper 
limits on the latent energy storage capacity, from either equation (22) or 
Figure 33, are 6923 kJ/m^ (610 Btu/ft^) and 17 114 kJ/m^ (1508 
Btu/ft^ ), respectively. 

Case b - Assume an imposed heat flux density of 4305.6 J/s-m* (1365.7 Btu/h-ft^ ). 


• Determine an acceptable geometry and the corresponding latent energy 
storage capacity. From either Figure 46 or 57, the fin height cannot 
exceed 1.27 cm (0.5 in.) for a half-fin spacing 0.127 cm (0.05 in.). From 
equation (22) or Figure 33, the corresponding latent energy storage 
capacity is 2247 kJ/m^ (198 Btu/ft^). 


Example 4 

Given; A phase change device is required to have a fin height of no more than 
5.08 cm (2.0 in.) and its base is to be exposed to a heat flux density of 2152.8 J/s-m^ 
(682.9 Btu/h-ft^ ); 


Case a — The designer wishes to use a fin thickness of 0.0127 cm (0.005 in.) to 
store as much energy as possible. 
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Determine the base temperature rise if the designer selects a half-fin 
spacing of 0.3175 cm (0.125 in). From Figure 48, the maximum base 
temperature rise is 



ATj, - 0.21 X 291.67 61.3 K (1 10.3°R) 


which is too large for most practical considerations. 

• Determine the base temperature rise for a half-fin spacing of 0.127 cm 
(0.05 in.). From Figure 48, the maximum base temperature rise is 


^ 0.095 X 291.67 ^ 27.7 K (49.9°R) 


Case b - Assume the designer selects a fin thickness of 0.0508 cm (0.02 in.) 

• Determine the base temperature rise for a half-fin spacing of 0.3175 cm 
(0.125 in.). From Figure 43, the maximum base temperature rise is 


ATf, == 0.088 X 291.67 ^ 25.7 K (46.2°R) 


• Determine the base temperature rise for a half-fin spacing of 0.127 cm 
(0.05 in.). From Figure 43, the maximum base temperature rise is 


AT^j 0.043 X 291.67 12.5 K (22.6°R) 


Case c — Assume a fin thickness of 0.127 cm (0.05 in) is selected. 

• Determine the base temperature rise for a half-fin spacing of 0.3175 cm 

(0.125 in.). From Figure 38, the maximum base temperature rise is 

AT 0.060 X 291.67 17.5 K (31.5°R) 


• Determine the base temperature rise for a half-fin spacing of 0.127 cm 
(0.05 in.). From Figure 38, the maximum base temperature rise is 


AT^ ^ 0.030 X 291.67 ^ 8.8K(15.8°R) 


Example 5 

Given: A phase change device is to be designed to have its base exposed to a heat 
flux density of 2152.8 J/s-m^ (682.9 Btu/h-ft^). Space constraints restrict the plan area 
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to 0.0557 (0.6 ft^) and the maximum base temperature at 5400 s (1.5 h) is not to 

exceed the PCM fusion temperature by more than 8.33 K (15°R). Consequently, 


8.33 

291.67 


28.6 X 10-3 


An estimate of the required latent energy storage capacity is 


E” = 2152.8 X 5400/1000 = 11 625 kJ/m^ (1024 Btu/ft^) 


Determine the acceptable geometry to satisfy these conditions. An examination of 
Figures 33 through 35 reveals that the maximum acceptable height is 6.35 cm (2.5 in.) 
and that the required height must be greater than this as the fin spacing is decreased and 
as the fin thickness is increased. 

Case a - Consider a fin thickness of 0.0127 cm (0.005 in.). The specified 
conditions cannot be met for this case (at least considering geometries used herein) as can 
be seen in Figure 54. 

Case b — Consider a fin thickness of 0.0508 (0.02 in.). Figure 55 shows that a fin 
height of 7-62 cm (3.0 in.) and a half-fin spacing of 0.127 cm (0.05 in.) will facilitate 
maintenance of the base temperature rise constraint. An examination of Figure 34 
indicates that the latent energy storage capacity is also satisfied. Figure 67 shows that the 
time required for the PCM to melt is 5760 s (1.6 h.). 


V. RELATED TOPICS 


In the preceding section, attention has been given to certain properties of several 
paraffins and to an extensive numerical study of a particular phase change device. The 
results of the parametric study should be helpful to the designer of phase change devices 
for applications compatible with the model, but the user needs an awareness of features 
that may require more specialized treatment in a critical design. Some of these important 
features which are discussed in this section include: (1) filler, (2) convection, (3) 
solid-phase thermal conductivity, and (4) consideration of an apphcation with a 
nonuniform thermal boundary condition. 


A. Filler 

The principal reason for using a filler is to improve the flow of energy to or from 
the PCM. Most nonmetallic PCM’s, particularly the paraffins, have very low thermal 
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conductivities; consequently, large thermal gradients may be required to etTect the desired 
heat transfer between the PCM and the heated or cooled surFaee of the phase change 
device. Incorporation of the metallic filler within tlie PCM region yields a composite 
medium with a higher effective thermal conductivity. If the application imposes a size 
limitation on the phase change device, utilization of a filler reduces the volume available 
for the PCM. Consequently, when considering the energy storage capacity, the large heat 
of fusion of the displaced PCM is replaced by the specific heat-temperature rise product 
for the filler. 

Fillers are discussed in References 2, 15, and 16. Some types which have been 
considered are: 


• Aluminum powder 

• Aluminum foam 

• Aluminum wool 

• Aluminum honeycomb 

• Copper foam 

• Alumina (AI 2 O 3 ) foam 

• Alumina (AI 2 O 3 ) powder 

• Aluminum fins (straight). 

Based on tests involving aluminum wool, aluminum foam, copper foam, and aluminum 
honeycomb, it is reported [15] that honeycomb offers distinct advantages over the other 
three types. The results of a number of performance tests on honeycomb are given. 

Straight fins offer some advantage over honeycomb in certain cases since better 
thermal contact between the filler and the housing of the phase-change device may be 
possible. Fins can be welded to the housing or they may be provided integrally with the 
housing via milling or casting. 

The numerical study treated in Sections III and IV applies specifically to 
straight-finned arrangements. Likewise the experimental work presented in Reference 2 
pertains to applications with straight fins. 

The designer should evaluate the advantages and penalties of using a filler. 
Fabrication, thermal contact resistance, added weight, reduced energy storage capacity, 
and improved effective thermal conductivity are all interrelated factors. The effects of 
straight-fin spacing, height, and width on the thermal performance of the model treated 
in Sections III and IV are reflected in the results. Some comparative insight into the 
effect of the use of honeycomb filler may be obtained from the results of tests presented 
in Reference 15. 
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B. Convection 


The physical model and the parametric numerical study which are presented in 
Sections III and IV are based on the conductive mode of heat transfer only within the 
PCM. If motion occurs within the liquid phase of the PCM, the heat transfer is affected. 
This behavior may require special consideration for applications where the liquid is 
subjected to surface or body forces of sufficient magnitude to initiate and sustain the 
motion. For the considerations given herein, there appears to be no previous or 
contemporary treatment of the formidable problem of an exact analysis of convective 
motion in enclosed cells of liquid. Some experimental correlations exist for certain cases, 
and these have been applied with some success to phase change problems. A brief 
background review and some considerations related to PCM studies are presented in the 
following paragraphs to aid the user in his study of convective effects. 

1. Brief Review of Convective Considerations. If one neglects body force fields 
other than gravity and surface forces other than surface tension, convective motion 
within the liquid of a PCM could still be due to one or more of the following effects: 

• Buoyancy forces 

• Surface tension forces 

• Density change as the PCM undergoes a phase change. 

In early studies of buoyancy driven convection, Lord Rayleigh [17] related the 
Nusselt numbers (Nu) to the product of the Grashoff number (Gr) and the Prandtl 
number (Pr). The product of the latter two is known as the Rayleigh number (Ra). 

Numerous experimenters have conducted studies that were related to 
hydrodynamic instability caused by buoyancy as it effects the heat transfer process, 
including the studies cited in References 18 through 29. O’Toole and Silveston [30] 
verified the Nusselt number versus Rayleigh number correlation for a fluid confined 
between two parallel plates (Fig. 108). A number of authors, including Edwards and 
Catton [31] later showed the effect of L/d ratios for closed cells (Fig. 109), where L is 
the cell height and d represents fin cell spacings. 

The critical Rayleigh number is defined as the value at which convection begins. 
In the region below the critical value, the heat transfer is substantially by conduction 
only (i.e., Nu= 1). The critical value, which depends on the boundary conditions, is 1708 
for a fluid confined between two infinite horizontal, isothermal, conducting walls; it is 
720 for the case of a fluid confined within a nonconducting walled container [29]. As 
shown for a bounded cell (Fig. 109a), this critical value tends to increase as the cell sides 
approach one another (i.e., as L/d increases), and it also tends to increase as the walls 
become more conducting (Fig. 109b). 

Although scientists have been aware of surface tension driven convection for some 
time, quantitative studies of this phenomena appear to be scarce. The flow patterns 
created in Benard’s classic experiment [32], which produced cellular circulation patterns 
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RAYLEIGH NUMBER. Ra 

Figure 108. Silveston’s experimental results in the neighborhood of instability 
in various liquids (Fig. 1 1 of Reference 2). 


b 1 

Figure 109. Region of influence of L/d on heat transfer (Fig. 12 of Reference 2). , 



INITIATION OF 
CONVECTION 


heat TRANSFER 
WITH 6% OF 
NO LATERAL 
WALL CASE 



t NATURAL CONVECTIC 

CwiTHIN NATURAL CONVECTION ; 

^5% OF CASE WITH >TXv DECREASED DUE TO 

if^suLATiNO lateral. 



1.0 2.0 

height-to-width ratio <L/d) 



149 



in a very shallow liquid, were initially attributed to buoyancy ert'ccts; however. Block 
[33], Pearson [34J, Scriven [35], and Sternling [36J later proved that this phenomenon 
was due to surface tension driven convection. The initial discovery of this plienomenon is 
attributed to Marangoni [37], and tlie term “Marangoni Flow” is commonly associated 
with this phenomenon. 

The nature of the Marangoni (low, as discussed by Young, Goldstein and Block 
[38], is that temperature variations across a free gas/liquid interface prt)duced variations 
in the shear force along the surface. This is due to the dependency of surface tension on 
temperature, which is estimated by Gambrill [39] to be linear; 


o - Oo + bT 


(23 I 


Since the coefticient “b” is negative, an increase in temperature at the surface is 
accompanied by a subsequent decrease in the surface tension. Hershey [40] has shown 
qualitatively that a depression occurs in the surface at a local hot spot, causing the liquid 
to flow away from the hot zone and toward the cold zone (Fig. 110). McGrcw and 
Larkin [41] photographed this effect for a number of configurations, producing dramatic 
verification of the phenomenon. 


HOT FLUID 




Figure 1 10. Surface tension convection patterns (Fig. 13 of Reference 2). 


The level of surface tension driven convection 
Marangoni number, Ma, which is given by 


is correlated by using the 


Ma = 


-da dT 
dT dy 


L2 


pna 


(24) 
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The critical Marangoni number for a fluid fixed between a rigid and free surface is given 
as 80. 


Nields [42] has stated that the onset of convection in a standard one-g field 
might be better determined by correlating the Nusselt number with a normalized 
parameter, R, given by 


Ra Ma 
-|- 

Ra^r ^^cr 


(25) 


This assumes that buoyancy and Marangoni driven convection are additive. 
Grodzka [43], however, noted that experimental data on buoyancy and Marangoni 
surface effects indicated that such an additive relation does not hold. Regardless, the 
convective currents caused by Marangoni flow in a standard one-g gravitational field are 
usually small compared to those caused by buoyancy. For most fluids at normal 
temperatures, Pearson [34] has shown that a liquid thickness of 1 cm or less must be 
attained before Marangoni effects overshadow buoyancy effects. 

For systems with small characteristic dimensions, the nondimensional Bond 
number. Bo, is given by 


Bo 


Gravity Forces _ pgl? 
Surface Tension Forces a 


(26) 


This number is sometimes used to evaluate the relative importance of Marangoni effects 
as compared to buoyancy effects. From equation (26), it follows that a low Bond 
number indicates a high degree of surface tension effects. 

Although surface tension effects exist at all unlike interfaces, only the liquid/gas 
interface (as opposed to liquid/liquid and liquid/solid interfaces) is expected to produce 
appreciable resulting flows. However, no proof of this conjecture was noted in the 
literature surveyed. 

Volume change driven convection can be caused by the phase change process. 
During freezing, the new layer of frozen material at the interface tends to contract, since 
the solid density is usually greater than the liquid density. Consequently, the liquid near 
the front will flow toward the interface to fill the volume shrinkage caused by 
solidification. During melting, the liquid at the interface tends to flow away to allow for 
the volume created by melting. Tien and Koump [44] have stated that this effect will 
cause both the freezing and melting process to be retarded. This retardation is due to the 
ingress of warm fluid during freezing and the egress of warm fluid during melting (Fig. 

Ill) 
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Figure 111. Volume change driven flow patterns (Fig. 14 of Reference 2). 

If liquid circulation is created in the melting processes, warm fluid could be 
drawn into the interface, thereby, augmenting the melting process. 


Tien and Koump [44] have also shown, in a computational exercise for a 
fictitious system where the solid density is 25 percent greater than the liquid density, 
that only a 10 percent reduction in the freezing rate occurs due to volume change 
effects. Since paraffins experience only a 5 to 10 percent volume increase on melting, it 
can be inferred that an even smaller effect can be expected in paraffins for similar 
conditions. 


In summary, the principal modes of heat transfer that are involved in a typical 
one-g static phase change process are conduction and convection. Convection is caused 
only by buoyancy driven currents for reasonable container sizes. Surface tension effects 
are negligible, except for very thin films, and volume change effects may be ignored. 

2. Zero Gravity . Since a number of applications of phase change devices [2] 
have been in space flight situations, some attention is devoted here to zero-g 
performance. 

When comparing environs of near Earth orbital space with that of the Earth, 
several differences are noted in the heat transfer process. In Earth orbit, the reduced 
gravitational force is nearly balanced by the centripetal orbital force, creating an effective 
zero-g environment. Effects caused by reduced pressure, radiation field, meteoroid 
bombardment, and three-dimensional spacecraft maneuvers are also possible. 


Restricting the hypothetical phase change device under consideration to be a 
hermetically sealed container that is isolated from exterior thermal effects by insulation 
and antipenetration shields and to be aboard a nonmaneuvering vehicle, then reduced 
gravity remains as the only important alien effect. 

The primary effect of reduced gravity on the heat transfer mechanism is in the 
lessening or elimination of buoyancy convection. Typical measurements of the net 
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gravitational acceleration force on a spacecraft indicates acceleration levels of the order 
1 X 10"’ g. Using nonadecane paraffin properties, a typical Rayleigh number of this low 
level of gravitational acceleration is 


Ra = 0.01 AT (27) 

Equation (27) indicates that, for reasonable container sizes (i.e., less than 15 cm (6 in.) 
cell depths), the temperature difference across the liquidus portion of the cell must be 
^ 444°C (830° F) to produce buoyancy driven convection. Since this temperature 
difference is well in excess of normal operating values, buoyancy stimulated convection 
may be considered negligible. 

At one-g, Marangoni flow or surface tension-driven convection is normally 
unimportant; however, at zero-g this is not necessarily true. As discussed earlier, liquid 
flow caused by surface tension at a liquid/vapor interface may occur. 

A group of experiments, performed during the mission of NASA’s Apollo- 14, 
revealed that the surface tension driven phenomena in zero-g is a reality and can produce 
significant convection [45,46]. Using data from these experiments, Grodzka [46] has 
plotted the relation existing between the Marangoni number and the ratio of effective 
thermal conductivity to actual conductivity (Fig. 112). These data show a rapid increase 
in the convective level at a Marangoni number slightly greater than 300. The temperature 
difference in the Krytox test liquid was only 2.5°C (4. 1°F) in this instance. 

Close examination of these data stimulates some questions. Applying these data to 
nonadecane paraffin contained in a cell with a 15.24 cm (6 in.) characteristic dimension 
at a Marangoni number of 300, a Kgff/K of 12 is predicted in Figure 1 12. This is a very 

high convective level. The temperature difference across a parallel plate system required 
to reach this convective level is only 0.001 1°C (0.002°F); however, for this same system, 
a Rayleigh number of 4.32 X 10® is predicted at one-g. From this, a of only 2.5 

is indicated, which is a convective level well below that predicted by the tentative surface 
tension data. This conflicts with earlier investigators, who reported that for this condition 
at one-g, Marangoni driven convection is unimportant. This could be explained, at least 
partially, by the fact that the Apollo 14 data were taken for a free surface, whereas the 
Figure 108 correlations were taken for top and bottom bounded cells. 

Using nonadecane paraffin properties, evaluated at the phase change temperature, 
the Marangoni number is given by 


Ma = 2374 ATL 


( 28 ) 


This indicates that for reasonably sized cells, a temperature difference of only 0.003 K 
(0.006° F) is necessary for the onset of Marangoni driven convection in a paraffin filled 
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device at zero-g. This fact, along with the high levels indicated in Figure 112, suggests 
that Marangoni convection can be appreciable at zero-g. 

Since the magnitude of volume change effects discussed earlier were low for 
freezing (where buoyancy driven convection is negligible) as well as for melting, it may 
be implied that the volumetric effects are negligible for pure conduction as well as 
convective processes. Consequently, the volumetric effects on the heat transfer process 
are negligible for zero-g operation. Heat transfer modes occurring in a phase change 
device at zero-g are then reduced to Marangoni-driven convection and pure conduction. 
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Figure 112. Heat transfer characteristics of B^ard cells (Fig. 15 of Reference 2). 
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A secondary effect of reduced gravity wliicli can significantly alter the heat 
transfer process is the ullage gas [)osition. In a typical rectangular cell containing ullage 
volume, the ullage gas may configure itself iti any of a number of possible modes. A 
number ol possible cell ullage locations are given by Reference 2, those of more 
importance are shown in Figure 113. When heating/cooling from the bottom, only 
configurations d and e (Fig. 113) would alter the normal process; however, when 
healing/cooling from the top or sides, configurations a, b, d, and e would all reduce the 
rate of lieat transfer due to the insulating effect of the ullage gas. Small bubbles 
occurring in the lupiid could induce convective currents. In a zero-g field, this motion 
could be caused by the Marangoui How phcnomoii causing the bubbles to migrate toward 
warm zones. Also, bubbles could be entrapped in the freezing solid (configuration d of 
lug. 113), thereby decreasing tlie apparent thermal conductivity of the PCM. Fortunately, 
paraffins have the property of being good surface wetters, which iji zero-g tends to force 
the ullage to form in the center of the cell (configuration c of Fig. 113); however, 
insufficient (luantitative data are available on these phenomena to determine the effect on 
the heat transfer process under given conditions. 



I'igurc 1 13. Zero-g ullage configurations (Pig. 16 of Reference 2). 


155 






Grodzka [4] has examined effects of the space environment on the microscopic 
processes. She concluded that complex coupling effects between phase change kinetics 
and various possible complex coupling effects between phase change kinetics and various 
possible modes of convective motion cannot be predicted accurately without actual flight 
data. She also concluded that the magnitude of magnetic and electric fields likely to be 
encountered in the Earth’s orbit are not expected to alter phase change behavior 
significantly from that observed on Earth. Finally, she stated that radiation fields 
encountered in Earth orbit are expected to have little effect, except perhaps in the case 
of. an organic PCM where long time exposures will result in the buildup of impurity. 
Although the later statement applies to a paraffin, no definitive information on this 
effect could be found in the literature surveyed. 

Although insufficient reduced gravity data are available to corroborate these 
flndings, some observations have been made as to the mechanisms of heat transfer in a 
phase change device. When a free surface is present and thermal conditions are proper, 
Marangoni convection may be present. Insufficient data are currently available to 
determine quantatively the influence of this phenomenon. In the absence of surface 
tension effects, conduction is the primary process governing heat transfer. The 
conduction process may, however, be altered by the ullage gas location. Furthermore, 
secondary effects of property degradation due to radiation may occur if the PCM is 
organic. 

V 3. PCM Related Considerations. Regarding buoyancy-driven convective currents 
in Ijie liquid during the melting of a PCM, a major factor is the positioning of the heated 
surface. If the heated surface is located at the bottom of a PCM cell during melting, 
convection may significantly augment the heat transfer. Conversely, the role of 
convection is minor, vanishing for one-dimensional cases, when the heated surface is at 
the top of the cell. The contribution of convection in the liquid to the heat transfer 
during solidification depends on positioning of the cooled surface, but it also depends on 
the amount by which the initial liquid temperature exceeds the fusion temperature of the 
PCM. 


Prensner and Hsu [47] reported on studies of natural convection in air within 
closed cells. Their work revealed the possibility of a number of different 
geometry-dependent convective patterns, even under steady state conditions. Such studies 
attest to the formidable task of analytically predicting the transport mechanism occurring 
during natural convection within a bounded fluid. When considering convection within 
the liquid during melting of a PCM, the problem is further complicated by the fact that 
one of the boundaries is moving (i.e., the interface). Even if accurate predictions of the 
convective mechanism were possible, another unpredictable feature occurring in finned 
PCM devices is the tendency for the PCM to melt up the fin with the solid PCM core, 
subsequently, falling to the heated surface. 

While the aforementioned difficulties appear pessimistic, some success has been 
reported [2, 48] on approximately accounting for convection during the melting of a 
paraffin. The approach made use of the correlations given by O’Toole and Silveston [30] 

which are : 
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1700 < Ra < 3500 

Nu = 

3500 < Ra < 10® 

Nu = 

10* < Ra < 19® 

Nu = 


0.00238 

0.229 Ra°-^ (29) 

0.104 (Ra)° -3°5(Pr)« ®«‘’ 


Although the correlations given by equation (29) were determined from data on steady 
state heat transfer through horizontal layers of fluid, they were used as a means of 
estimating the convective effect within the liquid of finned PCM cells similar to the one 
modeled in Section 111. The scheme involved continuous replacement of the thermal 
conductivity of the liquid by an effective thermal conductivity, which was calculated by 


Kgff = Nu X k 


(30) 


It was recognized that this scheme would not yield true temperatures in the liquid 
because the convective mechanism therein was being artificially replaced by conduction in 
a hypothetical medium with a larger thermal conductivity than that of the liquid PCM. 

Figure 1 14 shows some experimental data for the interface location during 
melting of a paraffin inside a cell similar to the one modeled in Section 111. Also shown 
are numerical predictions for conduction only and for inclusion of convection in tiie 
manner previously described. The numerical predictions shown in Figure 114 were 
achieved with an explicit formulation which was subject to a stability requirement. For 
nodes located in the fin, this requirement demanded excessive computer time; 
consequently, temperatures along the fin were approximated with the aid of 
measurements from thermocouples attached to the fins. With the approximate inclusion 
of convection, the predictions were promising. The predicted interfacial positions 
followed the same trend as the experimental data. Although magnitudes did not exactly 
agree, the difference appeared to be due to a shift in time and might have been 
associated with starting conditions. Use of this approach for including convection is 
discussed more extensively in Reference 2. Also presented is a more detailed discussion of 
the computer program and several comparisons between predictions and experimental 
data. The computer program is also listed and briefly discussed in Appendix D. 

4. Solid Phase Thermal Conductivity . Nearly all published thermal conductivity 
values for paraffins apply to the liquid phase. Since the mechanism of heat conduction is 
typically different in solids than in liquids, one should expect some difference between 
the thermal conductivity value for the solid phase and that for the liquid phase of a 

PCM. 


In studies reported in References 2 and 49, attempts to analytically predict the 
transient interfacial location during the phase change of a paraffin were more successful 
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FUSION FRONT LOCATION AT CELL 



Figure 1 14. Comparison between predicted and measured fusion front location 

for an example melt run. 
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for the case of melting than for the case of solidification. Two plausible explanations for 
the lack of success with solidification existed. One explanation was simply that the true 
value of the thermal conductivity of the solid phase may not have been used in the 
analysis. Knowledge of the thermal conductivity of the solid is more critical in analyses 
of solidification than in those of melting. Another possible explanation was related to the 
character of the solid-liquid interface during the phase transition. During melting, the 
interface appeared smooth. During solidification, however, it was covered with numerous 
projections called dendrites. Their presence may have affected the interfacial heat transfer 
rate. The discrepancy may also have been due to a combination of both effects. 

Thomas and Westwater [50] studied the nature of the interface during freezing 
and melting of n-octadecane. Measured interfacial velocities exceeded analytical 
predictions by as much as 100 percent; the difference was attributed to irregularities at 
the interface. 

Bailey and Davila [51] monitored the position of the solid-liquid interface during 
solidification of hexadecane and octadecane. Bailey and Liao [52] found that effective 
solid thermal conductivities, which yielded agreement between predictions and 
experimental measurements, exceeded published values by factors ranging from 2.5 to 3. 

Dyer and Griggs [49] measured interfacial positions during solidification of 
nonadecane and hexadecane; they also determined the values of solid thermal 
conductivity which forced agreement between the measurements and numerical 
predictions. Tne work was extended to include measurement of the thermal 
conductivities of solid nonadecane and hexadecane. The measured values did not agree 
with the values required to force agreement between experimental and analytical results 
(Table 16). 

As previously mentioned, values of thermal conductivity for the solid phase of 
paraffins are less prominent in the literature than are those for the liquid phase. Some of 
the available values are summarized in Table 16 for reference. 

In summary, analytical predictions of the solidification of paraffins suffer 
accuracy due to the aforementioned effects. The designer needs to exercise some 
judgement in the selection of an appropriate thermal conductivity value to use for the 
solid phase. Uncertainty associated with the effect of interfacial irregularities (i.e., 
dendrites) renders the problem less amenable to an exact treatment. 


C. Example Application with a Nonuniform Thermal 
Boundary Condition 

The parametric study, which is presented in Section IV, pertains to a numerical 
analysis of a symmetrical half of a single PCM cell. The results can be extended to 
multicell phase change devices by appropriate scaling provided that the thermal boundary 
condition imposed on the full-scale device is uniform and agrees with the type analyzed. 
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TABLE 16. SOME REPORTED THERMAL CONDUCTIVITY VALUES FOR THE 
SOLID PHASE OF CERTAIN PARAFFINS IN W/m-°C 
(Btu/h-ft-°F) 



a. Value used to force agreement between analytical and experimental results. 


One application of phase change devices involves heat transfer between the device 
and a flowing fluid. In such applications, the thermal condition existing at the separating 
surface between the PCM and the fluid is not necessarily uniform. An analysis of the heat 
transfer between the fluid and the PCM is not a simple problem. The entire apparatus 
including PCM, housing, and flow passage could be subdivided into a thermal network of 
nodes and interconnecting thermal resistances, and an appropriate numerical study could 
be performed. Such an approach, however, could involve a large number of nodes, 
particularly if much internal detail is needed. For a specific case, this approach may be 
the most desirable as well as productive; however, it certainly does not appear feasible to 
attempt to generalize it into an extensive parametric study, not only because of the vast 
nature of the task but also because of the large number of parameters which would be 
involved. The unwieldiness of the problem is evident from the integral analysis presented 
by Baily and Liao [55] even where some simplifying assumptions were incorporated. An 
approximate approach which utilizes input from the parametric study is outlined in the 
following paragraphs. 

The model under consideration is depicted in Figure 1 15. Fluid flows through the 
flow passage and makes intimate thermal contact with a surface separating it from a 
device filled with a PCM. For this model, it is assumed that the lower side of the flow 
passage is insulated or that a symmetrical PCM device is also located on the other side. 
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Figure 115. Physical model for heat transfer between a phase change 
device and a flowing fluid. 


An energy balance on a differential length of the flow passage yields the following 
differential equation: 


ax 


p A„c 


m 


p p at 


+ m c 


ax 

pT 


a^x 




(31) 


Neglecting any effect due to the surface separating the PCM and fluid, the instantaneous 
heat transfer between the surface and the PCM can be equated to the instantaneous heat 
transfer between the flowing fluid and the surface. Xhis equality can be expressed by 


h(Tm-Ts) = Uo(X^-Xf) 


(32) 


Xhe right-hand side of equation (32) represents the instantaneous heat transfer per unit 
area between the surface and the PCM. Xhe two-dimensional character of the problem has 
been masked by use of an overall heat transfer coefficient based on the difference 
between the surface temperature and the PCM fusion temperature. Equation (32) can be 
rearranged to yield the surface temperature as 
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(33) 


‘m 


f 


1 + (Uq/H) 1 + h/Uo 


Equation (33) implies that — Tj^ when h greatly exceeds Uq and Tg — Tf when Uq 
greatly exceeds h. 

The elimination of Tg in equation (31) with the aid of equation (33) and the 
neglection of heat conduction in the flowing fluid give 


dT 


m 


3T 


+ m c 


m 


P 3x 


= hB 


m 


^1+h/Uo 1+h/Uo 


) 


(34) 


Equation (34) can be solved numerically for both spatial and temporal variations of Tj^; 

the solution, however, requires input of the convective heat transfer coefficient h for the 
flowing fluid and the overall coefficient Uq for heat transfer with the PCM. The 

computer programs (listed in Section IV) provided a value of Uq as part of the output. 

This value is determined by dividing the instantaneous heat transfer per unit area between 
the surface and the PCM by the difference in the surface temperature and the fusion 
temperature of the PCM. Although the computations apply to a symmetrical cell, use of 
the Uq values may be acceptable since lateral heat conduction between adjacent PCM 

cells is probably much less significant that heat transfer between the surface and the 
PCM. If this hypothesis is valid, the Uq values may serve as a convenient aid in analyzing 

the variation of T^^. The solutions discussed in the following paragraphs are restricted to 
those cases for which the aforementioned assumptions are meaningful. 

The time derivative can be approximated by 

at At ^ 


The space derivative can be approximated by 

^^m ^ ^*’^i+l,n ^*^i+l,n-l 
dx Ax 


The first subscript denotes time while the second subscript denotes space. 


(36) 
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Insertion of the difference expressions given in equations (35) and (36) into 
equation (34) and some rearrangement result in 


^i+l,n 


1 


pCpA mcp 

^ + h/Ug) AT ^ At 

At Ax 1 + h/UQ 


T 


rn 


iqi 

(37) 


Equation (37) facilitates computation of the temperature of a node in the fluid at time 
t + At from the previous temperature of that node at time t and the temperature of the 
adjacent node at x - Ax and at time t + At. The nodal pattern is illustrated in Figure 1 16. 
The computations require input of the initial temperature distribution in the fluid (i.e., 
temperatures along the horizontal axis of Figure 116) and the transient variation of the 
fluid temperature at the inlet to the flow passage (i.e., temperatures along the vertical 
axis of Fig. 116). The internal convective heat transfer coefficient for the flowing fluid 
must be estimated from knowledge of the character of the flow (e.g., laminar, 
fully-developed, turbulent, etc.). The value of Uq is dependent on the amount of PCM 

which has undergone phase change; consequently, the computational scheme should 
incorporate sequential alteration of Uq. One method employed in the following 

computations involved representing Uq as a function of the fraction of PCM which has 

undergone phase change. The fraction is determined at each location and used in the 
functional representation to determine Uq for use in the computation of T^^. 



(x distance INDEX) 

Figure 1 16. Nodal pattern for calculating bulk temperature of flowing fluid. 
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A computer program which facilitates computations of the type previously 
outlined is listed in Appendix C. Exploratory calculations have been made, but the 
predictions have not been corroborated by any experimental data. The program was 
developed and intended as a tool for investigating the effect of the phase change process, 
flow rate, and internal convective heat transfer coefficient on the temperature at the exit 
of a flow passage which exchanges energy with the phase change device. Without some 
history of comparison with experimental data, the integrity of the results can only be 
conjectured at this stage. 

The following parameters were used in two computer runs to examine the 
influence of the internal convective heat transfer coefficient inside the flow passage. The 
results were generated for the case of a sinusoidal variation of the fluid temperature at 
the inlet of the flow passage. Initially, the flow was stagnant with the fluid in the passage 
and the entire phase change device considered to be at the fusion temperature of the 
PCM. 

• Length of phase change device — 31.7 cm (1.04 ft) 

• Width of phase change device — 7.6 cm (3 in.) 

• PCM thickness — 3.8 cm (1.5 in.) 

• Flowing fluid — Water 

• Flow rate ~ 1.89 gm/s (15 Ib/h) 

• PCM fusion temperature - 32.2°C (90'^F) 

• Uq variation with fraction of PCM melted, f - Uq = 86.02 - 76.45 f 
Btu/h-ft^-°F or Uq = 488.36 - 434.03 f J/s-m^-°C 

• Temperature variation of fluid at passage inlet - = 110 + 10 sin 

(6.283 t) °F with t (h) or Tj^iet = 43.33 + 5,56 sin (2.26 X 10^ t) °C 
with t (s) 

• Convective heat transfer coefficient inside flow passage — 1st Run: 
h = 283.86 J/s-m^-^C (50 Btu/h-ft^-^F); 2nd Run: h = 28,39 J/s-m^-^C (5 
Btu/h-ft^-°F). 

The results for the fluid temperature at the exit of the flow passage are shown in 
comparison to the inlet temperature variation in Figure 117. For both heat transfer 
coefficients, the temperature drop across the passage decreases as the PCM melts (as 
would be expected). The relative influence of the heat transfer coefficient on the 
temperature drop across the passage is evident. 



PRESCRIBED INLET FLUID TEMPERATURE 



do'3UniVU3dlAI31 
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Figure 1 17. Example exit temperature variations for a sinusoidal temperature at inlet. 






APPENDIX A 


PROPERTY APPENDIX 


Tabulations of certain properties for a number of phase change materials are given 
in this appendix. They have been included principally to serve as a convenient reference 
for the user; they are not intended to be exhaustive. For more detail regarding the 
advantages and disadvantages of the listed materials, the referenced literature should be 
consulted. 

Following is an index for the tabulations: 

Table A-1. Fusible Materials with a Heat of Fusion Greater than 185.96 kJ/kg (80 
Btu/lb) Listed in Order of Increasing Melt Temperature from 4.4 to 65°C 
(40 to 150°F). (Table 3-1 of Reference 15). 

Table A-2. Comparisons of Fusion Properties of Selected PCM’s Reported in Reference 
56 (From Table 3 of Reference 56) 

Table A-3. Properties of Six PCM’s Studied in Reference 57 with a Fusion Temperature 
Between 273 and 373 K (From Table 2 of Reference 57) 

Table A-4. Melt Temperatures of Fifteen Low Temperature PCM Candidates (From 
Table 6 of Reference 57) 

Table A-5. Four Low Temperature PCM’s Recommended in Reference 57 (From Table 
7 of Reference 57) 

Table A-6. PCM’s Used in Tests Reported in Reference 58 (From Table 1 of Reference 
58) 

Table A-7. Comparative Data on Waxes and Inorganic Hydrates (From Reference 59) 

Table A-8. Melting Temperature of Some Mixtures Investigated for Air Conditioning 
Thermal Energy Storage (From Table 3-1 of Reference 1) 

Table A-9. Melting Temperature and Heat of Fusion of Inorganic Salt Hydrates (From 
Table 3-3 of Reference 1) 

Table A- 10. Melting Point, Heat of Fusion, and Latent Heat Density of Some Inorganic 
Hydrates (From Table 3-4 of Reference 1) 

Table A- 11. Melting Temperature and Heat of Fusion of Some Inorganic Hydrate 
Eutectics (From Table 3-5 of Reference 1) 
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Table A-12. Calorimetric Data for Some Organic Thermal Energy Storage Materials for 
Solar Heating (From Table 3-6 of Reference 1) 

Table A-13. Melting Temperature and Heat of Fusion of Organic Materials for Air 
Conditioning-Hydrocarbons and Halogenated Hydrocarbons (From Table 3-7 
of Reference 1) 

Table A-14. Melting Point and Heat of Fusion for Some Organic Eutectics (From Table 
3-7 of Reference 1) 

Table A-15. Melting Temperature and Heat of Fusion for Some Acetamide-Organic 
Eutectics and Compounds (From Table 3-7 of Reference 1) 

Table A-16. PCM Data for Some Alcohols, Phenols, Aldehydes, Ketones, and Ethers 
(From Table 3-7 of Reference 1) 

Table A-17. PCM Data for Organic Acids and Miscellaneous Organics (From Table 3-7 of 
Reference 1) 

Table A-18. Calorimetric Data for Air Conditioning Organic Materials (From Table 3-8 of 
Reference 1) 

Table A- 19. Melting Point and Heat of Fusion for Clathrates and Semi-Clathrates Melting 
Congruently Above 0°C (32°F) (From Table 3-9 of Reference 1) 

Table A-20. Melting Point for Organic-Inorganic Eutectics (From Table 3-10 of Reference 

1 ). 
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TABLE A-1 . FUSIBLE MATERIALS WITH A HEAT OF FUSION GREATER THAN 
185.96 kJ/kg (80 Btu/lb) LISTED IN ORDER OF INCREASING MELT 
TEMPERATURE FROM 4.4 TO 65°C (40 TO 150°F). (FROM 
TABLE 3-1 OF REFERENCE 15). 


Material 

Melting Point 
°C (°F) 

Heat of Fusion 
kJ/kg (Btu/lb) 

Tetradecane 4 H 3 q 

5.6 (42) 

227.6 (98) 

Formic Acid HCOOH 

7.8 (46) 

246.1 (106) 

Pentadecane Cj 5 Hj 2 

10 (50) 

206.7 (89) 

Myristic Acid Ethyl Ester 
CH 3 (CH 2 ), 2 C 00 C 2 H 5 

11 (51) 

185.8 (80) 

Acetic Acid CH 3 CO 2 H 

17 (62) 

185.8 (80) 

Hexadecane C, ^ H 3 4 

18 (64) 

236.8(102) 

Lithium Chloride Ethanolate 
UCP 4 C 2 H 60 

21 (69) 

185.8 (80) 

n-Heptadecane Cj 7 H 3 ^ 

22 (71) 

213.6 (92) 

d-Lactic Acid CH 3 CHOHCOOH 

26 (79) 

185.8 (80) 

Octadecane C j ^ H 3 g 

28 (82) 

243.8 (105) 

13-Methyl Pentacosane C 26 H 54 

29 (84) 

195.1 (84) 

Methyl Palmitate C 2 7 H 3 4 O 2 

29 (84) 

104.3 ( 88 ) 

Nonadecane Cj 9 H 40 

32 (90) 

220.6(95) 

Trimyristin (Cj 3 H 2 7 COO )3 C 3 Hj 

33 (91) 

202.0 (87) 

2-Dimethyl-n-docosane C 2 4 H 5 0 

35 (95) 

197.4 (85) 

Eicosane C 20 H 42 

37 (98) 

246.1 (106) 

1-Tetradecanoi CH 3 (CH 2 )i 2 CHj OH 

38 (100) 

229.9 (99) 

Camphenilone C 9 Hj 4 O 

39 (102) 

204.3 ( 88 ) 

Caprylone (CH 3 (CH 2 )a )2 CO 

40 (104) 

257.7(111) 

Docosyl Bromide C 2 2 H 4 5 BR 

40 (104) 

202.0(87) 

Heneicosane C 2 1 H 44 

41 (105) 

213.6 (92) 

7-Heptadecamone C ^ 7 H 3 4 0 

41 (105) 

199.7 ( 86 ) 

1-Cyclohexyloctadecane C 2 4 H 2 s 

41 (106) 

218.3 (94) 

4-Heptadecanone C 17 H 34 O 

41 (106) 

197.4(85) 
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TABLE A-1. (Concluded) 


Material 

Melting Point 
°C (°F) 

Heat of Fusion 
kJ/kg (Btu/lb) 

8-Heptadecanone Cj 7H34O 

42 

(107) 

202.0 (87) 

Cyanamide CH2N2 

44 

(111) 

209.0 (90) 

Docosane C2 2 H4 ^ 

44 

(112) 

248.5 (107) 

Methyl Eicosanate C2 1 H4 2 O2 

45 

(113) 

227.6 (98) 

Tricosane C23H48 

47 

(117) 

232.2 (100) 

3 -Heptadecanone Cj 7 H3 4 0 

48 

(118) 

216.0 (93) 

2 -Eptadecanone C 1 7 H3 4 0 

48 

(119) 

216.0 (93) 

Camphene C^oHi ^ 

50 

(122) 

239.2(103) 

9 -Heptadecanone Cj 7H34O 

51 

(123) 

211.3 (91) 

Tetracosane C2 4 H5 0 

51 

(124) 

253.1 (109) 

Elaidic Acid Cj 3H34O2 

51 

(124) 

218.3 (94) 

Methyl Behenate C24H46O2 

52 

(126) 

232.2(100) 

Pentacosane C2 5 H5 2 

54 

(129) 

236.8(102) 

Ethyl Lignocerate C2 e H5 2 O2 

54 

(129) 

216.0(93) 

Hypo Phosphoric Acid H4P2O5 

55 

(131) 

213.6 (92) 

n-Hexacosane €75^54 

56 

(133) 

255.4(110) 

My ristic Acid C 1 3 H2 7 COOH 

57 

(135) 

199.7 (86) 

Heptacosane C2 7 H5 ^ 

59 

(138) 

234.5 (101) 

Ethyl Cerotate C2 gHs ^02 

60 

(140) 

222.9 (96) 

Octacosane €28^53 

61 

(142) 

253.1 (109) 

Nonacosane C2 9 0 

64 

(147) 

239.2(103) 

Stearic Acid 7 H3 5 CO2 H 

64 

(148) 

199.7 (86) 

Triacontane C3 0 Hg 2 

65 

(150) 

250.8 (108) 
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TABLE A-2. COMPARISONS OF FUSION PROPERTIES AND MELT TEMPERATURE OF 
SELECTED PCM’S REPORTED IN REFERENCE 56 (FROM TABLE 3 AND MELT 
TEMPERATURE OF REFERENCE 56) 
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.. As reported in Reference 56. 






















TABLE A-3. PROPERTIES OF SIX PCM’S STUDIED IN REFERENCE 57 WITH A FUSION TEMPERATURE 
BETWEEN 273 AND 373 K (FROM TABLE 2 OF REFERENCE 57) 
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TABLE A-4. MELT TEMPERATURES OF FIFTEEN LOW TEMPERATURE PCM 
CANDIDATES (FROM TABLE 6 OF REFERENCE 57) 


PCM 

Melt Temperature 
K (°R) 

n-Dodecane 

263.4 (474.1) 

Transit-Heat 

257.2 (463.0) 

Formic Acid/Formamide (75.5/24.5) 

256.4 (461.5) 

Transit-Heat 

252.2 (414.0) 

Acetic Acid/Formic Acid (53.1/46.9) 

250.7 (451.3) 

Acetic Acid/Water (58%) 

246.8 (444.2) 

n-Decane 

243.3 (437.9) 

Formic Acid /Potassium Formate (75/25) 

241.7 (435.1) 

Heptanone-4 

240.2 (432.4) 

Hydrazine Hydrate (N2H4 * H2O) 

233.2 (419.8) 

Monomethylhydrazine 

220.8 (397.4) 

Tricaproin 

213 (383.4) 

Tributyrin 

198 (356.4) 

Ammonia 

195 (351.0) 

Ammonia Hydrate (NH3 * H 20 j 

194 (349.2) 
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TABLE A-7. COMPARATIVE DATA ON WAXES AND INORGANIC HYDRATES (FROM REFERENCE 59) 
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TABLE A-8. MELTING TEMPERATURE OF SOME MIXTURES INVESTIGATED FOR 
AIR CONDITIONING THERMAL ENERGY STORAGE (FROM 
TABLE 3-1 OF REFERENCE 1) 


Composition 

Mole 

Weight 

(gram) 

Relative 

Weight 

(%) 

Melting 
Point 
°C (°F) 

NA2SO4 

142.04 

37.2 


10H2O 

180.15 

47.4 

18(64) 

NaCl 

58.44 

15.4 


Total 

380.63 

100.0 


Naj SO4 

142.04 

37.8 


10 H2O 

180.15 

48.0 

10(50) 

NH4CI 

53.49 

14.2 


Total 

375.68 

100.0 


Naj SO4 

142.04 

37.6 


IOH2O 

180.15 

47.6 

13 (55) 

1/2 NaCl 

29.22 

7.7 

I/2NH4CI 

26.75 

7.1 


Total 

378.16 

100.0 


Naj SO4 

142.04 

35.8 


IOH2O 

180.15 

45.4 

7(45) 

KCl 

74.55 

18.8 


Total 

396.74 

100.0 
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TABLE A-9. MELTING TEMPERATURE AND HEAT OF FUSION OF INORGANIC SALT 
HYDRATES (FROM TABLE 3-3 OF REFERENCE 1) 
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TABLE A-9. (Concluded) 
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TABLE A-10. MELTING POINT, HEAT OF FUSION, AND LATENT HEAT DENSITY 
OF SOME INORGANIC HYDRATES (FROM TABLE 3-4 
OF REFERENCE I) 


Material 

Melting Point 
°C (°F) 

Heat of Fusion 
(Literature) 
kJ /Btu\ 

Latent Heat Density 
kJ /Btu\ 

m^VTP'/ 

CaClj • 6 H 2 O 

28.9 (84) 

170.2 (73.3) 

2.86 X 10* (7688) 

MgClj • 6 H 4 O 

117.2 (243) 

172.3 (74.2) 

2.69 X 10* (7226) 

Ca(N03)2 ■ 4 H 2 O 

42.8 (109) 

142.1 (61.2) 

2.59 X 10* (6954) 

Mg(N03>2 6 H 2 O 

95.0 (203) 

159.8 ( 68 . 8 ) 

2.34 X 10* (6271) 

Zn(N 03>2 • 6 H 2 O 

36.1 (97) 

130.0 (56.0) 

2.69 X 10* (7219) 

A1(N03)3 9 H 2 O 

70.1 (158) 

155.6 (67)^ 



a. Experimental value 
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TABLE A- 11. MELTING TEMPERATURE AND HEAT OF FUSION OF SOME 
INORGANIC HYDRATE EUTECTICS (FROM TABLE 3-5 
OF REFERENCE I) 


Composition 

Wt.% 

Melting Point 
°C (°F) 

Heat of Fusion 
kJ /Btu \ 

) 

CaCIj 

37 



Ca(N 03>2 

23 approx. 

13 (56) 


HjO 

40 



CaClj 

41 



MgCl2 

10 

25 (57) 

174.5 (75.1) Experimental 

H2O 

49 



CaCIj 

48.0 



KCI 

4.3 

27 (81) 


NaCl 

0.4 



HjO 

47.3 



Ca(N 03>2 -4H2 0 

45 

25 (77) 

129.7 (55.8) Calculated 

2 n(N 03)2 • 6 H 2 O 

55 



Ca(N 03>2 • 4H2O 

67 

30 ( 86 ) 

136.0 (58.5) Calculated 

Mg(N 03>2 - 6 H 20 

33 



Ca(N 03>2 -4H2 0 

72 

35 (95) 

138.9(59.8) Calculated 

AI(N03)2 ■ 9 H 2 O 

28 



Mg(N 03)2 • 6 H 2 O 

53 

61 (142) 

148.1 (63.8) Calculated 

A1(N03)2 • 9 H 2 O 

47 



Mg(N03)2 • 6 H 2 O 

18 

32 (90) 

130.1 (56.0) Calculated 

2 n(N 03)2 • 6 H 2 O 

82 




I8I 






TABLE A- 12. CALORIMETRIC DATA FOR SOME ORGANIC THERMAL ENERGY STORAGE MATERIALS 

FOR SOLAR HEATING (FROM TABLE 3-6 OF REFERENCE 1) 
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TABLE A- 13. MELTING TEMPERATURE AND HEAT OF FUSION OF ORGANIC 
MATERIALS FOR AIR CONDITIONING-HYDROCARBONS AND 
HALOGENATED HYDROCARBONS (FROM TABLE 3-7 
OF REFERENCE 1) 


Compound 

Melting Point 
°C (°F) 

Heat of Fusion 
kJ/Btu\ 

i^viry 

n-Tetradecane 

5.9 (42.6) 

225.5 (97.1) 

n-Pentadecane 

9.0 (48.2) 

162.9 (70.1) 

n-Hexadecane 

18.2(64.8) 

233.9(100.7) 

n-Hexadecylbromide 

14(57.2) 

— 

1-Heptadecane 

11.2(52.2) 

131.1 (56.5) 

1-Octadecane 

17.9(64.2) 

129.1 (55.6) 

1-Pentadecyne 

10(50) 

— 

1-Hexodecyne 

15 (59) 

■ — 

p-Xylene 

13.2(55.8) 

159.9 (68.9) 

Tridecyibenzene 

10(50) 

168.3 (72.5) 

Tetradecylbenzene 

16 (60.8) 

173.3 (74.6) 

2-Nonyl naphthalene 

12(53.6) 

— 

3-Chlorobiophenyl 

17(62.6) 

— 
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TABLE A-14. MELTING POINT AND HEAT OF FUSION FOR SOME ORGANIC 
EUTECTICS (FROM TABLE 3-7 OF REFERENCE 1) 


Compound 

Mole % 

Melting Point 
“C (°F) 

Heat of Fusion 
kJ /Btu\ 
kg \ lb / 

a Chloroacetic Acid 
Phenol 

31.0 

69.0 

16.5(61.7) 

144.1 (62.1) 

/3 Chloroacetic Acid 
Phenol 

35.0 

65.0 

11.8(53.2) 

138.2(59.5) 

a Chloroacetic Acid 
p Cresol 

28.0 

72.0 

13.4(56.1) 

141.6 (61.0) 

a Chloroacetic Acid 
o Cresol 

30.9 

69.1 

15.8(60.4) 

189.6 (81.7) 

Acetic Acid 
Benzoic Acid 

88.0 

12.0 

9.5 (49.1) 

184.6 (79.5) 

Naphthalene 
Diphenyl Methane 

25.0 

75.0 

14.5 (58.1) 

109.0 (46.9) 
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TABLE A-16. PCM DATA FOR SOME ALCOHOLS, PHENOLS, ALDEHYDES, 
KETONES, AND ETHERS (FROM TABLE 3-7 OF REFERENCE 1) 


Compound 

Melting Point 
°C (°F) 

Heat of Fusion 
kJ / Btu\ 

kg VIT/ 

Alcohols and Phenols 



1-Decanol 

6.1 (43.0) 


l-Undecanol 

15 (59.0) 


Glycerol 

18 (64.4) 

200.4 (86.2) 

m-Cresol 

12 (53.6) 

94.0 (40.5) 

2-Chlorophenol 

7 (44.6) 


3,3 Dimethyl Cyclohexanol 

11 (51.8) 


Aldehydes, Ketones, and Ethers 



1 ,4 Dioxane 

11.8(53.2) 

145.3 (62.6) 

2-Decanone 

14 (57.2) 


n-Tridecanol 

15 (59.0) 


Parahlehyde 

12.6 (54.7) 

104.4 (45.0) 

2,5 Dichloroaceiophenone 

14 (57.2) 
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TABLE A-17. PCM DATA FOR ORGANIC ACIDS AND MISCELLANEOUS 
ORGANICS (FROM TABLE 3-7 OF REFERENCE 1) 


Compound 

Melting Point 
°C (°F) 


Organic Acids 



Formic Acid 

8.4 (47.1) 

276.7 (1 19.0) 

Acetic Acid 

16.6 (61.9) 

201.3 (86.6) 

Acrylic Acid 

12.3 (54.1) 

154.9 (66.6) 

n-Octanoic Acid 

16.3 (61.3) 

148.2(63.7) 

n-Nonanoic Acid 

12.3 (54.1) 

128.1 (55.1) 

Oleic Acid 

14 (57.2) 


Other 



Diethyl Isophthalate 

1 1.5 (52.7) 


Dimethyl Adipate 

8.5 (47.3) 


1,2 Ethylene Diamine 

8.5 (47.3) 


m-Cresyl Acetate 

12 (53.6) 

_ 

2,5 Dimethyl Analine 

14.2(57.7) 
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TABLE A-18. CALORIMETRIC DATA FOR AIR CONDITIONING ORGANIC MATERIALS 

(FROM TABLE 3-8 OF REFERENCE 1) 
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TABLE A-19. MELTING POINT AND HEAT OF FUSION FOR CLATHRATES AND 
SEMI-CLATHRATES MELTING CONGRUENTLY ABOVE 0°C (32°F) 

(FROM TABLE 3-9 OF REFERENCE I) 


Composition 


Type I Clathrate Hydrates 
HjS • 6.1 H^O 
HjS-e.lHjO 
CO 2 6.0 H 2 O 
CI 2 • 7.3 H 2 O 
SO 2 • 6.0 H 2 O 
SO 2 • 6.1 HjO 
C 2 H 4 O • 6.9 H 2 O 

Type II Clathrate Hydrates 
C 4 H 8 O ■ 17.2 H 2 O 
(Tetrahydrofuran) 

C 4 H 8 O • 2 H 2 S ■ 17 H 2 O 
(CHj) 3 COH • 2H2S • 17H2O 

Amine Semi-Clathrate Hydrate 
(CH 3 ) 3 N • 10 1/4 H 2 O 

Tetrabutylammonium Salt Semi-Clathrate 

Hydrates 

BU 4 NF • 32 H 2 O 
BU4NCI • 32H2O 
Bu4NBr • 32 H 2 O 
BU4NNO3 • 32H2O 
BU4NOH • 32H2O 
BU4NHCO3 ■ 32H2O 
(Bu4N)2HP04 • 64H2O 
(Bu4N)2C2 04 • 64 H 2 O 
BU4NCHO2 • 32H2O 
BU4NCH3CO2 • 32H2O 
BU4NC6HSCO2 • 32H2O 

Tetraisoamylammonium Salt Semi-Clathrate 

Hydrates 

i-Am4NF -40 H 2 O 
i-Am 4 NCl • 38 H 2 O 
i-Ani 4 NOH • 40 H 2 O 
i-Am4NCH02 • 40 H 2 O 


Melting Point 
°C(°F) 



0.4 (33) 
29.5 (85) 

9.9 (50) 
9.7(50) 
7.0(45) 

12.1 (54) 

11.1 (52) 


246.9 (106.3) Calculated 


4.4(40) 255.2 ( 109.0) Experimental 

21.3 (70) 276.1 (1 18.9) Experimental 

7.3 (45) ) 


5.9 (43) 238.5 (102.7) Experimental 


24.9 (77) 

15.7 (60) 

12.5 (54) 

5.8(42) 

30.2(86) 

17.8(64) 

17.2(63) 

16.8(62) 

12.5 (54) 184.1 (79.3) Experimental 

15.1 (59) 209.2 (90. 1) Experimental 

3.5 (38) 


31.2(88) 
29.8 (86) 

31 (88) 

15-20(59-68) 
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TABLE A-20. MELTING POINT FOR ORGANIC-INORGANIC EUTECTICS 
(FROM TABLE 3-10 OF REFERENCE I) 


Urea Eiiteetics-Binary 








Melting Point 

A 



Mole A 



Wt. % A 


°C(°F) 

NH 4 CI 



17.5 



15.9 


101 (214) 

NIl4Br 



23,5 



33.4 


76(169) 

NH 4 1 



25.5 



45.2 


66(151) 

NH 4 NO 

3 


47.5 



54.7 


46(115) 

NH 4 SCN 


39.5 



48.4 


32 (90) 

LiNOj 



16 



18 


76(169) 

LiNOj 



45.5 



48.1 


104 (219) 

NaNOa 



23 



29.5 


84 (183) 

KNO 3 



16 



24.3 


110(230) 

CaCNOa 

)2 


9.4 



22 


92(198) 

Ca(N03 }2 


25 



48 


140t286) 

Sr(NOa 

>2 


15 



38 


67 (153) 

Ba(NOa 

h 


10.5 



33.8 


96 (205) 

KaCOa 



7 



14.8 


102 (216) 

Urea Eutectics-Ternary 









Melting Point 

A 

B 

Mole % A 

Wt. 

77 A 

Mole % B 

Wt. % B 

°C (°F) 

NH 4 CI 

NH4NOa 

4 


2.4 

46 

53 

35 (95) 

NH 4 C! 

NH 4 SCN 

4 

2.5 

40 

46 

20 (68) 

NH 4 NO 3 

LiNOa 

34 

40 

7 

7.1 

34 (93) 

NH 4 NO 3 

LiNOa 

46 

51 

32 

31 

65 (149) 

NH 4 NOa 

NaNO 

3 

41 

47 

8 

9.7 

35 (95) 

NH4NOa 

KNOa 


44 

49.9 

4 

5.7 

43 (109) 

LiNOa 

NaNO 

3 

24 

26 

10 

13 

56 (133) 

L 1 NO 3 

NaNO 

3 

45 

49 

3.5 

4.6 

100(2121 

LiNOa 

KNOa 


16 

17 

6 

9.5 

62(144) 

LiNOa 

KNOa 


45 

43 

19 

27 

96(205) 

LiNOa 

Ca(NOa )a 

16.5 

18.1 

1.5 

4 

62 (144) 

LiNOa 

Ca(NOa )a 

41 

44.5 

4 

9.2 

82 (180) 

LiNOa 

Sr(NO 

3)2 

12 

11 

8 

23 

45 (113) 

LiNOa 

'lalNOa )a 

16 

17 

1 

4 

68(154) 

LiNOa 

Ba(NOa )a 

44 

46 

1 

4 

92 (198) 

NaNOa 

KNOa 


22.3 

27.7 

7 

10.3 

73 (163) 

NaNOa 

Ba(NOa >2 

17.5 

20.5 

4 

14.4 

68 (154) 

KNO 3 

Ca(NOa >2 

3 

4 

14 

20 

89(192) 

KNOa 

Ca(NOa >2 

18 

20 

24 

42.5 

125 (257) 

KNO 3 

BalNOj )2 

8.5 

1 1.5 

5.5 

19.3 

80 (176) 
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TABLE A-20. (Concluded) 


Acetamide Eutectics 




Melting Point 

A 

Mole % A 

Wt. % A 

°C(°F) 

NH4NO3 

32 

39 

38 (100) 

LiNOa 

20 

23 

25 (77) 

NaNOj 

15 

20 

59 (138) 

KNO 3 

6 

i 

9.8 

72(162) 
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APPENDIX B 


THERMAL RESISTANCE VARIATION 


Some of the results of the parametric study which are shown in Section V exhibit 
interesting trends. Specifically, some plots of the temperature rise at the termination of 
melting versus fin height demonstrate a relative maximum, a relative minimum, or both. 
Several curves in Figures 5 1 mrough 62 indicate this behavior. A plausible explanation of 
this behavior is presented in this appendix. 

The physical model, which was treated in Section IV, is shown in Figure B-1. 
Heat flows directly from the bottom into the PCM. Heat also flows into the PCM 
through the fin and top. The relative thermal resistances of these paths is dependent on 
the geometry. 



Figure B-1. Physical model. 
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Some insight into the thermal resistance behavior may be gained by means of a 
simple but nonexact viewpoint. Suppose the effect of the top is ignored and the PCM 
and fin are each represented by a single node. The resulting thermal network, including 
thermal resistances, for the heat flow between the bottom and the PCM is indicated in 
Figure B-2. The equivalent thermal resistance for this arrangement is related to the 
parameters by 


KgqK^B= 1.0/[2W/H + 2 (K/K^)HS/(2H' +WS(K/K^)] (B-1) 


The right-hand side of equation (B-1) was evaluated for the range of geometrical 
parameters covered in the parametric study. The thermal conductivity ratio (K/K^) of 

1069 was also used. The results are plotted in Figures B-3 through B-5. The effect of 
geometry on the equivalent thermal resistance for the simplified network is shown by the 
curves. 


BOTTOM NODE 


H/KS2B W/2K^HB 

FIN NODE 


PCM NODE 


VW 

H/2K^WB 

Figure B-2. Simplified thermal network. 
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EQUIVALENT THERMAL RESISTANCE x K 



0 2 4 6 8 

PCM HEIGHT (in.) 


Figure B-3. Variation of equivalent thermal resistance with height for a 
fin thickness of 0.0127 cm (0.005 in.). 
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EQUIVALENT THERMAL RESISTANCE x K 













APPENDIX C 


COMPUTER PROGRAM FOR EXAMPLE DESCRIBED IN SECTION V 


An example application of a phase change device with a nonuniform thermal 
boundary condition is discussed in Section V. The equation used to compute the 
temperature variation of the flowing fluid is equation (37). The computer program shown 
in Table C-1 was written to facilitate the calculation of the mean fluid temperature from 
equation (37). This program occupies approximately 7 K words of storage. The data 
shown in Table C-2 is required for each run. In addition, the user must supply routines in 
the function subroutines COEFF and TEMP. COEFF calculates Uq and TEMP generates 

the inlet temperature of the fluid. Note that there are provisions for up to three different 
functions for Uq and the inlet temperature. The function used during any run is chosen 

by the data entered in the data card (Table C-2). There are provisions for multiple runs 
within one computer job. The program reads nine cards, processes the run, and then 
returns to read nine more. If none are present, the program terminates. For each run, 
header information is printed, such as flow rates, and fluid and PCM properties. Two 
time steps are printed simultaneously, one on the left-hand side of the page, the other on 
the right-hand side of the page. A description of the output data is given in Table C-3. 
At the end of each run, the total energy stored, QTOT, is printed. 
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TABLE C-i. COMPUTER PROGRAM USED TO COMPUTE ELUID TEMPERATURE 
VARIATION AS IT FLOWS THROUGH A PHASE CHANGE DEVICE 


real LENUTMi fiLUMP* moot 

INTEGER lin£#TlMECT 

OImENSIUN I MtANl 2i500) #QTAB(50’0) lOwOai 500»< TLUMPI 2 
1#500> 

COHMON 1BUI‘F/LINCT#0L0TME/N0DECT 
9998 READ ( 5/ lOU^ tND»9999 I TFREEZ 
100 FORMAT 

READ (SillU) SPHEAT#M00T»HTXFER# HI^T^»UtNBTH/XSTEF 
1#TIMELN# TMtsrp 

110 FORMAT ‘aFlO.2/2Fi0i2*F10»5^Fiat2<Fl0»5l 
read (5#16U) RHOiXSECT 
150 FORMAT '2F10*2| 

READ (5#19U) SPAC# VISC^ PRAN#COMD 
190 FORMAT (4F1U.2) 

READ (5#16U) MEIQHT#HTFU8/6PMLMP#C0NLi 1P*RH0LMP 
160 FORMAT ‘ M 0 . b# 4F1 0 . 2 » 

READ (5il7U) IrRIT*18UPX 
170 FORMAT (2H0) 

read (5/18UI ITMPINjIUORTN 
180 FORMAT (2U0; 

WRITE (6/19S(( 

199 format ( ' 1 ■ » 

WRITE l«*#2U0) 

200 FORMAT l ’ 0 ’/ T49/ • 

l¥¥^^¥(Mf#¥¥¥¥.v if im tf 4. ^ l^ tf tf if ifif. tf. t j 

WRITE <6/220( LENGTHjXSECT 

220 FORMAT (' 1 38/ ' CAPACITOR ♦ LENQTH ••F/»3/' FT*'»T8C 

l/'FLOw HASb.AREA -'F7.4/» 80.FT. ' , TllS/ ' * * » 

WRITE (6/230) HEIGHT/SPAC 

230 FORMAT '■ '/ 138/ ' DIMENSIONS ♦ PCM HEI3HI -'/F7.3/t 

1 FT* ‘/TBO, 'FLOWPASS HEIGHT ■*/F7«4,* »ir • ' / T115/ ' * ' ) 
WRITE (6/2'tOi WIDTH 

240 FORMAT (' '/149,'¥ H-WlOTH •'F7#4,' FT • • / TH5/ * ♦ • I 

WRITE (6/210) 

210 format '' '/ T 49/ ' ¥¥♦¥¥¥¥¥•¥♦♦¥¥¥¥¥♦»¥*♦**¥♦*♦¥♦¥*¥¥¥¥¥ 
!¥¥¥¥¥¥¥¥¥¥¥«¥¥¥¥¥¥¥¥¥¥¥¥¥¥¥¥¥¥*) 

WRITE (6*200) 

WRITE (6/2»0l MDOT/VISC 

250 FORMAT (' '/T38/'FLUID * PlOW RATE -'/F7.3/' LB 

1/HR*'/T8U ^'viscosity "'/FTiA/' l8FFT*HR* ' / T115/ ’ ¥ • ) 
WRITE (6/260) SPHEAT/PRAN 

260 FORMAT (' '/ 1 38/ ' PROPERTIES ¥ flPECIFji'HT -'/F5.3/' 

1 BTU/LB-I- ' /TbO/ *PRAND NO ■ ' / F7 • 3/ Tl I S* ' ♦ ' > 

WRITE (6*2^0) RHO/CONO 

270 FORMAT (' '/T49/'¥ DENSITY -•/F7»3/* .BFCUB FT'/T80 

l/'THERM LOND -'/F7.3/' BTU/HR-FT-F • / T115* ' * ' ) 

WRITE (6/210) 

WRITE (6/200) 

WRITE (6/280) HTFUS/CONLMP 

280 FORMAT (' ’/138/'PCM ¥ HEAT OF' FJS •'/F5*0/‘ 

1 BTU/CUB F I ' / T80/ * THERM CONO ■•/«'7*3/' 8 TU/MR-F T-F:' 



TABLE C-1. (Continued) 


2>TU5# '*' ) 

WRITE (6>2S#0) SPHLMP#TFREEZ 

290 FORMAT l* T 38* » PROPERTIES * iPECIFIS' MT -'/FSta*' 

1 BTU/LB-^ ' /TBO* 'FREEZE TEMP Fl' * T115# ' ♦ » ) 

WRITE (6*300) RHOLMP 

300 FORMAT (' ’* TA9* ' DENSITY 1.3'Coa FT'/T115*' 

1 ** ) 

WRITE (6*210) 

WRITE (6*200) 

WRITE (6*310) XSTEP*ITMPIN 

310 FORMAT (' '*T38*'RUN * K-STEP •'#F/»4** FlT»*T8C 

1* 'TEMp-iN MTN ! '*I1*T115* ’♦» I 
WRITE (6*320) TMESTP* lUORTN 

320 FORMAT (' ’* T38* ' PARAMETERS ♦ TIME-ST£P -'*F7.3*' MRS 
1*T80* *UO"R>N ! ' I1*T115» '*' » 

WRITE (6*330) TIMELN 

330 FORMAT '*TA9*'¥ TIME LENGTH ■* F7 . 3* » MRS ’ * T115* ' ♦ * 
WRITE (6*210) 

CI*RH0<*XSECt»8PHEAT/TMESTP 

C2»M00T»8PH£aT/X8TEP 

NODECT-LtNOTM/XSTEPAZ 

timect-timlln/tmestp 

VOLUME-wioIH'»XSTEPpMEIQHT 
MLUMP-R'^OLMp# VOLUME 
UC0NST-1/(MEI(3HT/( 2¥C0NLMP) ) 

0IA-2>¥SMAL 

REYNO-OIA¥M[)OT/(RHO¥XSECT¥VIsC» 

IHT8W*1 

IF (HTXFfcR«EQ.0«0) IHT8W-0 
DO I J-1*NUD£CT 
TMEAN( 1*0 >"TFREEZ 
TLUMPd* J)-TFREEZ 
I QTAB(JI-o* 

IRITSW-IWRIT-1 

IBUFF-0 

LINCT-25 

0LDTME-339»0 

00 3 TIMfe.-l*TlMECT 

T-TIMEPIMESTP 

IF ( IRITbW'EU* IWRIT) IRXTSW-0 

IRlTSw-IMI (Sw+l 

TMEAN(2*l )-TtMPlN( 1TMPIN*T) 

00 4 N00t-«2*N0DECT 
X-(NOOE-l I^XSTEP 

UO-COEFF IQI AB( NODE) *VOLUME*MTFUB*FRAC‘rN*IUORTN) 

IF (FRACTN'LT. 0.991 QO TO 3C 

TLUMP < 2* NOOE ) »OLOQ ( node » / < MLUMP*SPHLM»i> ♦TLUMP (1 * NODE ) 
TMPPCM*TLUMP(2 *n00EJ 
UO^UCONP I 
QO TO 40 

30 TMPPCM-TFRtEZ 
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TABLE C-1. (Continued) 


40 IP (lHT8M*tQ*i) 30 TO 20 
QRATZ-X/ » 01A*REYN04PRAN ) 

HTXPER- ( tONO/OI A ) 4HEAT < QRATZ > 

20 C4»(HTXFtR*wlDTH*U0»/(UO4HTXFER) 

CS"1/(C1+C2+C4» 

TMEAN(2iNOOe)-C5*(C44TMPPCM+C2*TiEAN( 2#NOWE-l)+ca 
1»TmEANU«NUDE) ) 

QAiiC2*( TnEAN(2«N00EI*TrlEAN< 2 «nOO£»1) I 
aB*Cl«(TnEAN(2«NOOE)oTMEANU«NOOEO ) 
Q—14(QA+QB)»ThE8TP*X8TEP 
QTAB( NODE >-QTa8( node >4Q 

0LDQ(N0Ot)-Q 

IP ( lRIT8M*eQ*tMRlT) CALL PRnTI T ( T*1£A Y>( 2# NOOE ) « QTAB 
1IN0DE)|FKACTN« J0#TlTRCAN(2«l)4Xf S00c4'iTXPibR#I6UPX) 

4 CONTiNUfc 

00 2 J-l«NUOECT 
TLUMP(1#J)-TLUMPI2#J) 

2 rMEAN(l#J)-Tf1EANl2#J> 

3 continue 
OTOT-O 

DO 10 J-2/N0DECT 
iC QTOT*QTOIt>UTAB( J) 

WRITE «6<140) QTOT 
140 FORMAT t&X#'WTOT *'»FlO»2) 

00 TO 98BP 
9999 CALL EX11 
END 

subroutine PRnTITITHEAN#QTA8#F^CTs» J 3»TiTEMPlN#X#N0DE 
l#HTXPR* I8UPXJ 

dimension T MEANl ( SOO ) « QTABI ( 500 ) « FRCTNl < 300 ) «U01 ( 500 ) 
1#HTXFRU800) 

common aBUFF^LINCT/OLOTME/NODECT 

IF ( ISUPX*tO.O.ANO»NOOE*NE»NOOECT» RETURN 

IP (PRCTN.OT'1.0) FRCTN»l*0 

IP (lBUFF»tQtO) 80 TO 10 

IF (OLOTME'EU.T) 30 TO 20 

IBUPP-0 

OLOTME-i 

00 TO 30 

20 TMEANl(NOOt)"TMEAN 
QTAB1(NU0E)-QTAB 
FRCTNUNUOt )-FRCTN 
U01(N00E)«UC 
HTXPRUNODti-HTXFR 
RETURN 

10 IF (0L0TMC*Ey«T) 00 TO 40 
IBUPP*! 

OLOTME-T 

prnttm-uloime 

PRNTTP-TtMPiN 
00 TO 20 
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TABLE C-l. (Continued) 


30 IF (LINCI iUE*6S) 30 TO 50 

WRITE (6MU0) PRNTTM>PRNTTP#T#TE>1RIS 
100 FORMAT I 'O' /T9, tTIME • ' # F6» 3, T88» ' TEM»I ••#F7.2#T84 
1#’TIME -'if^P'ajTlOai'TEMPlN ii'*F7.2) 

WRITE {6*\^0l 
00 TO 40 

50 WRITE I6>H0) PRNTTM#PRNTTP*T#TEMPIW 
no FORMAT ( '1 ' /T9, 'TIME ■ ' i F6» 3# T38» • T£M»iIN *«#F7.2>TB4 
1#'TIME #»‘6*3*T103i 'TEMPIN ■'«F7.2) 

WRITE (6*140) 

140 FORMAT ( ' 0 ' i I 9, • X • , Tl6* * TEMP ♦ , T28# ' 3 • » T36, • CHAN3E ' * T47 
1* 'U0ST56* 'HTXFER’ #T80# 'X'*T88* 'TEMP*»T100, '(iSTl09 
2* 'CHANGE' * T119i 'U0'*Tl26 * 'HTxFER' I 
LlNCT-4 

40 IF (LINCT 30 TO 60 
WRITE 16*160) 

160 FORMAT I ' I ' * T 9* ' X ' * T1 6* ' TEMP ' # T28* • 3 ' * T36* ' CHAN3E ' * T47 
1* 'UO' /T55< 'HTXFER' *T80* 'X'*T88* ' TEMP' »T10W» • Q' # T109 
2* 'CHANGE'* 1119* 'UO'* Tl 26* 'HTXFER' ) 

LINCT-2 

60 LINCT-HNCI+1 

WRITE (6*160) X*TMEANHN00E)*QTa31(N03C»*''RCTN1(N00E) 
1*U01(NOOEJ *hT XFRl ( NODE I * X* TmEAN* QT A3 » FRCTN* UO/ HTXFR 
150 FORMAT (' ' * 6Fl 0 • 2/ 11 X* 6FlO • 2 ) 

RETURN 

END 

function CUEFF ( QTAB* VOLUME* HTFUSiFR acts* I UORTN) 
VOLMLT-wi AB/HTFUS 
FRACTN*AB8( VOLMLT/VOLUME) 

IF (FRACTN»Q£*0*99) RETURN 
GO TO (1*2*3)*IU0RTN 

1 continue 

♦ U8EK SUPPLIED UO VARIATION ROUTINE' 

*** GOES BETWEEN STATEMENTS IS** 

**¥ selected when IUORTN-1* 

IF <FRACrN*QT.,25» GO TO 10 
COEFF-lub. 12-242. 566 FRACTN 
RETURN 

10 C0EFF-54*32-35.37*FRACTN 
RETURN 

2 CONTINUE 

*¥» A second UO VARIATION ROUTINE 
♦*¥ MAY GO BETWEEN STATEMENTS 263* 

*♦4 8ELECTEI3 when IUORTN-2* 

C0EFF-5UU.U 

RETURN 

3 CONTINUE 

¥»» A TMIRO UO VARIATION ROUTINE CAN 30 
*** AFTER statement 3* 

*** selected when IUORTN-3* 

RETURN 
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TABLE C-l. (Concluded) 


END 

function TEmPINI ITMPIN#T) 

QO TO (a*2/3l#ITMPlN 
I CONTINUE 

*** USER SUPPLIED INLET TEMPERATJRE 
*** ROUTINE GOES BETWEEN STATEMENTS 
♦#¥ 152* selected when XTEMPIN-I* 

IF (T«QT»l»5t 30 TO 10 
TEMPIN-1‘38*T 

RETURN 

10 IF IT.GT'S'SJ 30 TO 20 

TEMPIN-*‘^*8'»^*12»927WT+**826*T*T 

RETURN 

20 IF {T»a"lU*5l 30 TO 30 

TEMPIN*80S'0*S1N( (3»1*H5/16>»T> 

return 

30 IF <T.ar»lA.O) 30 TO 40 

TEMP IN-14SS. 165-1 57 •533pT^3»8SS*T»T 

RETURN 

40 TEMPIN*PI4«7-40*98*T 

return 

2 CONTINUE 

C *** A SELONO temperature VARIATION 

C *** ROUIANE MAY GO BETWEEN STATEMENTS 

C *** 253* selected when ITEMPIN-2‘ 

TEMPiN-IUO'O+SOiOvT 
RETURN 

3 continue 

c ¥*¥ A THiRu Inlet temperature variatun 

c *** ROUTINE MAY GO AFTER STATEMENT 

C 3* selected when lTEMPlN«3f 

RETURN 
END 

FUNCTION HEATcsRATZ) 

IF <aRATi*aT» iOl 1 30 TO lO 
EX— 0*26 

HEAT-2 •‘i/»QRATZs*EX 

RETURN 

10 IF (GRAT^.GT* .031 30 TO 20 
ex— 0«0J 

HEAT-7 •2E-3RATZ-SEX 
RETURN 

20 HEaT-7»54i 

RETURN 

END 
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1 ABLE C-2. REQUIRED INPUT DATA FOR THE PROGRAM LISTED IN TABLE C-I 
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TABLE C-2. (Concluded) 
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TABLE C-3. DESCRIPTION OF OUTPUT DATA FROM PROGRAM 

LISTED IN TABLE C-I 


Data Name 

Description 

Units 

TEMP IN 

Fluid inlet temperature 

°F 

X 

X location in the capacitor 

ft 

TEMP 

Temperature at X 

°F 

Q 

Accumulated heat transfer into the PCM 
at X 

Btu 

CHANGE 

Fraction undergone change of phase 


UO 

Uo 

Btu 

h-ft^-°F 

HTXFER 

Convective heat transfer coefficient between 
fluid and flow passage wall 

Btu 

h-ft^-°F 
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APPENDIX D 


COMPUTER PROGRAM USED TO FACILITATE THE EXPLICIT NUMERICAL STUDY 


The computer program used to facilitate the explicit numerical study is outlined 
and described in this appendix. A copy of the program for a melting run with inclusion 
of convective effects is given in Section D.l; the corresponding notation is described in 
Section D,2; and the steps are discussed in Section D.3. 

After listing this program, it was discovered that the program described in Section 
D.l would not run for the case of M = 1. The statements causing this incompatibility are 
indicated in Section D.l by an arrow placed at the left of the appropriate statement. 

A copy of the computer program for a solidification case is given in Section D.4. 
A comparison of the appropriate statements with tliose shown in Section D.3 generally 
illustrates the modification to accommodate M = 1. 

The statements wliich were included in the solidification program to allow 
for changing the effective interfacial area for nodes experiencing melting are underlined. 
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D.l. COMPUTER PROGRAM FOR MELTING 


1 

2 

3 

4 

5 

6 

7 

8 

9 

10 
11 
12 

13 

14 

15 

16 

17 

18 

19 

20 
21 
22 

23 

24 

25 

26 

27 

28 

29 

30 

31 

32 

33 

34 

35 

36 

37 

38 

39 

40 

41 

42 

43 

44 

45 

46 

47 

48 

49 

50 

51 

52 
5 ? 
54 


DIMENSION 

DIMENSION 

DIMENSION 

DIMENSION 

COMPUTATIONAL 

N=27 

M=4 

AM=M 

MM=M+1 

NN=N+2 


RV(5,31) , RH(5,31) , C(5,31) , QS(5,31), 
T2(5,31) , F(5,31) , Q(5,31) , QRAT(5.31). 
TM1(9), TM2(9), TM3(9), TM4(9), TIM(9) 
QRA2(5, 31) 

PARAMETERS 


Tl(5,31) 
T3(l ,31) 


NJ=N+3 

NI=N+4 

ND=9 

NDP=8 

TAU=1 .1 

KC0UNT=1 

MC0UNT=1 

MFIN=500 

EPS=1.E-06 

KCHK=1000 


J0E=1 


DT=1 .OE-04 

TAU2=(10.*DT)+(DT/3.) 
PHYSICAL PROPERTIES 
TAMB=80 . 

TIN=73.5 

HT0P=l.E-08 

HB0T=5.0 

6=(32. 2*3600. *3600.) 
WAX 

DEN=47.2 


TK=0.087 

CP=0.5 

TTR=73.04 

HTR=22.108 

TMELT=89.8 

HMELT=73.357 

TREF=50.0 

BETA=0. 00045 

VIS=14.3 


CL=CP 

TKL=TK 

BOTTOM PLATE DENOTED BY 1 
DEN1=171.0 
TK1=93.0 
CP1=0.22 

FIN DENOTED BY 2 
DEN2=171 .0 
TK2=93.0 
CP2=0.22 

TOP PLATE DENOTED BY 3 
DEN3=72.5 
TK3=0.09 


CP3=0.33 


207 



55 

56 

57 

58 

59 

60 
61 
62 

63 

64 

65 

66 

67 

68 

69 

70 

71 

72 

73 

74 

75 

76 

77 

78 

79 

80 

81 - 
82 

83 

84 

85 - 

86 

87 

88 

89 

90 ■ 

91 

92 

93 - 

94 

95 

96 

97 

98 

99 

100 
101 
102 

103 

104 

105 

106 

107 

108 
109 

no 


GEOMETRY PARAMETERS 


W=0. 75/12.0 


H=2. 625/12.0 
8 = 5 . 0 / 12.0 
S=W/(2.0*AM) 
$1=0.032/12.0 
$2=0.008/12.0 
$3=0.25/12.0 
F2JM=1. 0/(32.0*5*12.0) 


VERTICAL RESISTANCES 

RVO ,2)=(Sl/(TKl*S2*B))+(2./(HB0T*S2*B)) 
DO 10 1=2, MM 

10 RV(I,2)=(S1/(2.*TK1*S*B))+(1./(HB0T*S*B)) 
RV(1 ,3) = (S/(TK2*S2*B))+(S1/(TK1*S2*B)) 

DO 20 1=2, MM 

20 RV(I,3)=(1./(2.*TK*B))+(S1/(2.*TK1*S*B)) 
DO 30 J=4,NN 

30 RV(1,J)=((2.*S)/(TK2*S2*B)) 

DO 40 J=4,NN 
DO 40 1=2, MM 
40 RV(I,J)=(1./(TK*B)) 

RV( 1 ,N+3)=( S3/(TK3*S2*B) )+(S/(TK2*S2*B) ) 


DO 50 I =2, MM 

50 RV(I,N+3)=(S3/(2.*TK3*S*B))+(1./(2.*TK*B)) 

RV(1 ,N+4) = (S3/(TK3*S2*B))+(2./(HT0P*S2*B)) 

RV(2,N+4) = (S3/(2.*TK3*B*(S+(S2/2. ))))+(!. /(HT0P*B*(S+(S2/2, 


-^DO 60 I=3,MM 

60 RV(I,N+4)=(S3/(2.*TK3*S*B))+(1./(HT0P*S*B)) 
HORIZONTAL RESISTANCES 

RH(2,N+3)=((S2+S)/(2.*TK3*S3*B)) 

-►DO 70 I =3, MM 
70 RH(I,N+3)=(S/(TK3*S3*B)) 

DO 80 J=3,NN 

80 RH(2,J)=(S2/(2.*TK2*S*B))+(1./(2.*TK*B)) 


DO 90 J=3,NN 
—►DO 90 1 = 3, MM 
90 RH(I,J)=(1./(TK*B)) 

RH(2,2)=((S2+S)/(2.*TK1*S1*B)) 

-►DO 100 1=3, MM 

100 RH(I,2)={S/(TK1*S1*B)) 

RH(2,N+3)=RH(2,N+3)+RV(l ,N+3) 

NODAL CAPACITANCES 

C(1 ,2) = ((DENl*Sl*S2*B*CPl)/2.) 

DO no 1 = 2 , MM 

no C(I,2) = (DEN1*S1*S*B*CP1) 

DO 120 J=3,NN 

120 CO ,J) = ((DEN2*S2*S*B*CP2)/2.) 

C(1 ,N+3) = ((DEN3*S2*S3*B*CP3)/2.) 

C(2 ,N+3) = ( S3*B*(S+(S2/2 . ) )*DEN3*CP3) 
—►DO 130 1=3, MM 
130 C(I ,N+3) = (DEN3*S3*S*B*CP3) 

DO 140 J=3,NN 
DO 140 1=2, MM 

140 C(I ,J)=(DEN*(S**2)*B*CP) 

-►Q1 = (C(3,4)*(TTR-TREF)) 

Q2=Q1+((DEN*(S**2)*B*HTR)) 


)) 


) 
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1 1 1 ^Q3=Q2+( ( C(3 ,4) )*(TMELT-TTR) ) 

112 Q4=Q3+((DEN*(S**2)*B*HMELT)j 

113 INITIALIZATION OF PERTINENT QUANTITIES 

114 TIME=0.0 

115 QWAX=0.0 

116 QBW=0.0 

117 QSW=0.0 

118 QTW=0.0 

119 QT0P=0.0 

120 QFTR=0.0 

121 QTTR=0.0 

122 V1=0.0 

123 V2=0.0 

124 DM0=0.0 

125 DO 150 J=2,NJ 

126 DO 150 1=1, MM 

127 150 T1(I,J)=TIN 

128 ►IF(TIN.LT.TTR)QSIN=(C(3,4)*(TIN-TREF)) 

129 ^IF(TIN.GT.TTR.AND.TIN.LT.TMELT) QSIN=Q2+((C(3,4))*(TIN-TTR)) 

130 ^ IF(TIN.GT.TMELT)QSIN=Q4+((C(3,4))*(TIN-TMELT)) 

131 DO 160 1=1, MM 

132 160 T1(I,N+4)=TAMB 

133 DO 170 J=3,NN 

134 DO 170 1=2, MM 

135 170 QS(I,J)=QSIN 

136 DO 180 1=1, MM 

137 Q(I,2)=0.0 

138 180 Q(I,N+3)=0.0 

139 DO 190 0=3, NN 

140 190 Q(1 ,J)=0.0 

141 DO 199 0=2 ,NI 

142 DO 199 1=1, MM 

143 F(I,0)=0.0 

144 199 QRAT(I,0)=0.0 

145 READ(5,11) (TMl(I), I=1,ND) 

146 READ(5,11) (TM2(I), I=1,ND) 

147 READ(5,11) (TM3(I), 1 = 1 ,ND) 

148 READ(5,11) (TM4(I), 1=1 ,ND) 

149 READ(5,11) (TIM(I), 1 = 1 ,ND 

150 11 F0RMAT(8F10.0) 

151 WRITE(6,22) TIME,W,H,N,M 

152 22 F0RMAT(1X,5HTIME=,E15.8,10X,2HW=,E15.8,10X,2HH=,E15.8,10X,2HN=,12, 

153 15X,2HM=,12) 

154 WRITE(6,33)Q1 ,Q2,Q3,Q4 

155 33 F0RMAT(1X,3HQ1=,E15.8,10X,3HQ2=,E15.8.10X,3HQ3=,E15.8,10X,3HQ4=,E1 

156 5.8) 

157 WRITE(6,44) 

158 44 F0RMAT(2X,4HI 0 ,5X,19HVERTICAL RESISTANCE ,6X,21HH0RIZ0NTAL RESIST 

159 ANCE,6X,17HN0DAL CAPACITANCE, 6X .11HTEMPERATURE,5X.7HQS( 1 ,0) ) 

160 DO 200 0=2 ,NI 

161 DO 200 1=1 ,MM 

162 IF(0.EQ.N+4) GO TO 1 

163 GO TO 2 

164 1 QS(1,0)=0.0 

165 RH(I,0)=1.E08 

166 C(I,0)=0.0 
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167 2 IF(I.EQ.1)RH(I,J)=1.E08 

168 IF(J.EQ.2.0R.J.EQ.N+3)QS(I,J)=0.0 

169 IF(I.EQ.1)QS(I,J)=0.0 

170 WRITE(6 ,55)1 ,J ,RV( I ,J) ,RH( I ,J) ,C(I ,J) ,T1 (I ,J) ,QS( I ,J) 

171 55 F0RMAT(1X,I2,1X.I2,5X,E15.8.10X,E15.8.10X,E15.8,6X,E16.8,4X,E15.8) 

172 200 CONTINUE 

173 DO 889 J=3,NN 

174 889 T3(1,J)=T1(1,J) 

175 COMPUTATION SECTION COMPUTATION SECTION COMPUTATION SECTION 

176 3 TIME=TIME+DT 

177 DO 210 J=3,NJ 

178 QRAT(1 ,J) = ((T1{2,J)-TI(1 ,J)VRH(2.J))+((T1(1 ,J-1)-T1(1 ,J))/RV(1 ,J) 

179 +((T1(1,J+1)-T1(1,J))/RV(1,J+1)) 

180 QRAT(MM,J)=((T1(M,J)-T1(MM,J))/RH(MM,J))+({T1(MM,J-1)-T1(MM,J))/RV 

181 (MM,J))+((T1{MM,J+1)-T1(MM,J))/RV(MM,J+1)) 

182 210 CONTINUE 

183 DO 211 J=3,NJ 

184 ^DO 211 1=2, M 

185 QRAT(I,J) = ((T1(I-1 ,J)-T1(I,J))/RH(I,J))+((T1(I+1 ,J)-T1(I,J))/RH(I+ 

186 1,J))+((T1(I,J-1)-T1(I,J))/RV(I,J))+((T1(I,J+1)-T1(I,J))/RV(I,J+1) 

187 ) 

188 211 CONTINUE 

189 212 DO 220 J=3,NJ 

190 DO 220 1=1, MM 

191 220 QS(I,J)=QS(I,J)+(QRAT(I,J)*DT) 

192 DO 240 1=2, MM 

193 240 T2(I,N+3)=Tl(I,N+3)+((QRAT(I,N+3)*DT)/C(I,N+3)) 

194 T2(l,N+3)=T2(2,N+3) 

195 DO 250 J=3,NN 

196 DO 250 1=2, MM 

197 IF(QS(I,J).LT.Q1)T2(I,J)=TREF+(QS(I,J)/C(I,J)) 

198 IF(QS(I,J).GE.Q1.AND.QS(I,J).LE.Q2) T2(I,J)=nR 

199 IF(QS(I,J).GT.Q2.AND.QS(I,J).LT.Q3) T2(I ,J)=nR+((QS( I ,J)-Q2)/C( I , 

J)) 

200 IF(QS(I,J).GE.Q3.AND.QS(I,J).LE.Q4) T2(I ,J)=TMELT 

201 WHEN GOING FROM MELT TO FREEZE OR VICE-VERSA CHANGE THE FOLLOWING CARD 

202 IF(QS( I ,J) .GT.Q3.AND.QS( I ,J) .LT.Q4) F(I ,J)=(QS(I ,J)-Q3)/(DEN*(S**2 

203 )*B*HMELT) 

204 IF(QS(I,J)GT.Q4) T2(I,J)=TMELT+((QS(I,J)-Q4)/C(I,J)) 

205 WHEN GOING FROM MELT TO FREEZE OR VICE-VERSA CHANGE THE FOLLOWING CARD 

206 IF(QS(I,J).GE.Q4)F(I,J)=1.0 

207 IF(QS(I,J).LE.Q3) F{I,J)=0.0 

208 250 CONTINUE 

209 SPECIFICATION AND/OR DETERMINATION OF FIN TEMPERATURES 

210 THE FOLLOWING DO LOOP ASSUMES FIN TEMPERATURES FOR ITERATION 

211 DO 255 J=3,NN 

212 255 T2(l ,J)=T3(1 ,J) 

213 DO 260 L=1,NDP 

214 IF(TIME.GE.TIM(L). AND. TIME. LE.TIM(L+D) GO TO 4 

215 260 CONTINUE 

216 4 FAC=(TIME-TIM(L))/(TIM(L+1)-TIM(D) 

217 T2(l ,2)=TM1(L)+( TM1(L+1)-TM1(L))*FAC) 

218 T2(l ,8)=TM2(L)+((TM2(L+1)-TM2(L))*FAC) 

219 T2(l ,15)=TM3(L)+((TM3(L+1)-TM3(L))*FAC) 

220 T2(l ,21)=TM4(L)+((TM4(L+1)-TM4(L))*FAC) 

221 DO 270 1=1 , MM 
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222 270 T2(I,2)=T2(1 ,2) 

223 UNSPECIFIED FIN TEMPERATURES DETERMINED BY STEADY STATE EQUATIONS 

224 256 MC0UNT=MC0UNT+1 

225 DO 280 J=3,NN 

226 T3(l ,J)=T2(1 ,J) 

227 IF(J.EQ.8.0R.J.EQ.I5) GO TO 5 

228 IF(J.EQ.21) GO TO 5 

229 T2(l ,J) = ((T2(1 ,J-1)/RV(1 ,J) ) + (T2{2 ,J)/RH(2 ,J) )+(T2(1 ,J+1 )/RV( 1 ,J+I 

230 )))/({1./RV(1,J))+(1./RH(2.JT)+(1./RV(1,J+1))) 

231 5 CONTINUE 

232 280 CONTINUE 

233 IF(MCOUNT.GT.MFIN) GO TO 8 

234 DO 281 J=3,NN 

235 DIF=T2(1 ,J)-T3(1 ,J) 

236 IF(ABS(DIF).GT.EPS) GO TO 256 

237 281 CONTINUE 

238 IF(J0E.EQ.2) GO TO 285 

239 J0E=J0E+1 

240 DO 888 J=3,NJ 

241 DO 888 I=r,MM 

242 888 QS(I,J)=QS(I,J)-(QRAT(I,J)*DT) 

243 DO 282 J=3,NJ 

244 QRA2(1 ,J) = ((T2(2,J)-T2(1 ,J))/RH(2,J))+((T2(1 ,J-1)-T2(1 ,J))/RV(1 ,J) 

245 ) + ((T2(l ,J+1)-T2(1 ,J))/RV(1.J+D) 

246 QRA2(MM,J)=((T2(M,J)-T2(MM,J))/RH(MM,J))+((T2(MM,J-1)-T2(MM,J))/RV 

247 (MM,J)1 + ((T2(MM,J+1)-T2(MM,J))/RV(MM,J+D) 

248 282 CONTINUE 

249 DO 283 J=3,NJ 

250 ^DO 283 1=2 ,M 

251 QRA2(I.J)=((T2(I-1,J)-T2{I,J))/RH(I,J))+((T2(I+1,J)-T2(I,J))/RH(I4 

252 1 ,J)) + ((T2(I,J-1)-T2(I.J))/RV(I,J))+{(T2(I,J+i;-T2(I,J))/RV(I,J+1) 

253 ) 

254 283 CONTINUE 

255 DO 284 J=3,NJ 

256 D0284I=1,MM 

257 284 QRAT(I,J)=(QRAT(I,J)+QRA2(I,J))/2.0 

258 GO TO 212 

259 285 MC0UNT=1 

260 DO 286 1 = 1 ,MM 

261 286 QRAT(I,2)=({T2(I,2)-T1(I,2))*C(I.2)/DT) 

262 DO 287 1=1, MM 

263 287 QS{I,2)=QS(I,2)+(C(I.2)*(T2{I,2)-T1(I,2))) 

264 DO 290 1=2, MM 

265 QBW=(T2(I,2)-T2(I,3))/RV(I,3)+QBW 

266 290 QTW=(T2(I,N+3)-T2(I,N+2))/RV(l,N+3)+QTW 

267 QBT=QBW+((T2(1 ,2)-T2(l ,3))/RV(l ,3)) 

268 DO 300 J=3,NN 

269 300 QSW=QSW+((T2(1 ,J)-T2(2,J))/RH(2,J)) 

270 DO 310 1=1 ,MM 

271 310 QT0P=((T2(I,N+3)-TAMB)/RV(I,N+4))+QT0P 

272 DO 311 J=3,NN 

273 311 QFTR=QFTR+QRAT(1 ,J) 

274 DO 312 1 = 1, MM 

275 312 QTTR=QTTR+QRAT(1 ,NJ) 

276 QWAX=(QBW+QTW+QSW)*DT+(QWAX) 

277 RATIO=QSW/QBW 
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278 


ERR0R= ( ( QBT - ( QBW+QSW+QTW+QTOP+QFTR+QTTR) ) *1 00 . 0 ) /QBT 

279 


DO 930 J=3,NN 

280 


DO 930 1=2, MM 

281 


V2=V2+(F(I,J)*(S**2)*B) 

282 

930 

CONTINUE 

283 


DAVG=((2.0*V2)/(W*B))*12.0 

284 


QMELT= ( ( V2- VI ) *DEN*HMELT ) /DT 

285 


V1=V2 

* 286 


V2=0.0 

287 


DELT=ABS{T2(1 ,2)-TMELT) 

288 


PR=(VIS*CL)/TKL 

289 


RA= ( ( DEN**2 ) *G*CL*BETA*DELT*( DAVG**3 ) ) / ( V IS*TKL*1 728 . 0 ) 

290 


IF(T2(1 ,2).LE.TMELT) RA=0.0 

291 


IF(RA.GT.1.E05)TK=(TKL*0.104*(RA**0.305)*(PR**0.084)) 

292 


IF(RA.GE. 3500,0. AND. RA.LE.1.E05) TK=(TKL*0.229*(RA**0.252)) 

293 


I F ( RA . GE . 1 700 . 0 . AND . RA . LT . 3500 .0 ) TK= (TKL*0 .00238*( RA**0 .816 ) 

294 


TKR=TK/TKL 

295 


DO 945 1=2, MM 

296 


IF(F(I,3) .LT.0.25) GO TO 945 

297 


RV(I,3)=(1 •/(2.*TK*B))+(S1/(2.*TK1*S*B)) 

298 

945 

CONTINUE 

299 


DO 946 J=4,NN 

300 


DO 946 1=2, MM 

301 


IF(F(I,J).LT.EPS) GO TO 946 

302 


RV(I,J)=(1./(TK*B)) 

303 

946 

CONTINUE 

304 


DO 947 1=2, MM 

305 


IF(F(I,NN).LT.0.75) GO TO 947 

306 


RV(I,N+3)=(S3/(2.*TK3*S*B))+(1./(2.*TK*B)) 

307 

947 

CONTINUE 

308 


DO 950 J=3,NN 

309 


IF(F(2,J). LT.0.25) GO TO 949 

310 


RH(2,J)=(S2/(2.*TK2*S*B))+(1./(2.*TK*B)) 

311 

949 

CONTINUE 

312 - 


DO 950 1=3, MM 

313 


IF(F(I,J).LT.EPS) GO TO 950 

314 


RH(I,J)=(1./(TK*B)) 

315 

950 

CONTINUE 

316 


IF(KCOUNT.EQ.KCHK) GO TO 6 

317 


KC0UNT=KC0UNT+1 

318 


GO TO 7 

319 

6 

WRITE(6,66) TIME, QWAX, ERROR 

320 

66 

F0RMAT(1X,5HTIME=,E15.8,10X,5HQWAX=.E15.8,10X,6HERR0R=,E15.8) 

321 


DO 313 J=3,NN 

322 


AJ=J-3 

323 


BJ=J-2 

324 


IF{F(2,J) .GT.F2JM) DFL0=(BJ*S*12.0) 

325 


DFIN=(AJ*S*12.0) 

326 


IF(F(2,J).LT.EPS) GO TO 314 

327 

313 

CONTINUE 

328 

314 

DO 315 J=3,NN 

329 


AJ=J-3 

330 


IF(F(MM,J).LT.1.0) DMID=((AJ*S)+(F(MM,J)*S))*12.0 

331 


IF(F(MM,J) .LT.l .0) GO TO 316 

332 

315 

CONTINUE 

333 

316 

WRITE(6,67) DFIN, DFLO, DAVG, DMID 
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334 67 FORMAT( 1 X ,5HDFIN= ,E1 5 . 8 ,10X ,5HDFL0= ,E1 5 .8,10X ,5HDAVG= ,E1 5 .8 ,10X ,5H 

335 DMID=,E15.8) 

336 DDOT=(DMID-DMO)/DT 

337 DMO=DMID 

338 WRITE (6, 68) RA, TKR, QMELT, DDOT 

339 68 F0RMAT(1X,3HRA=,E15.8,10X,4HTKR=,E15.8,10X,6HQMELT=,E15.8J0X,5HDD 

340 0T=,E15.8) 

341 WRITE(6,77) QBW, QSW, QTW, QBT, RATIO 

342 77 FORMAT{ 1 X ,4HQBW= ,E1 5 .8 ,3X ,4HQSW= ,E15 .8 ,3X ,4HQTW= ,E1 5 .8 ,3X ,4HQBT= ,E 

343 15.8,3X,6HRATI0=,E15.8) 

344 WRITE(6,88) 

345 88 F0RMAT(2X,1HI,2X,1HJ,10X,11HTEMPERATURE,10X,15HFRACTI0N MELTED, lOX 

346 .IIHENERGY RATE . lOX .13HENERGY STORED) 

347 DO 320 J=2,NJ 

348 DO 320 1=1 ,MM 

349 WRITE(6 ,99) I ,J ,T2( I ,J ) ,F( I ,J) ,QRAT( I ,J) ,QS( I ,J) 

350 99 F0RMAT(1X,I2,1X,I2,8X,E15.8,8X,E15.8,8X.E15.8,8X,E15.8) 

351 320 CONTINUE 

352 KC0UNT=1 

353 7 QBW=0.0 

354 QSW=0.0 

355 QTW=0.0 

356 QT0P=0.0 

357 QTTR=0.0 

358 QFTR=0.0 

359 DO 330 0=2 ,NJ 

360 DO 330 1=1 ,MM 

361 330 T1(I,J)=T2(I,J) 

362 J0E=1 

363 IF(TIME.LT.TAU) GO TO 3 

364 8 WRITE (6,111) MCOUNT 

365 111 F0RMAT(1X,I3) 

366 STOP 

367 END 
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D.2. DESCRIPTION OF PROGRAM NOTATION 


AJ 

AM 

B 

BETA 

BJ 

C(U) 

CL 

CP 

CPI 

CP2 

CPS 

DAVG 

DDOT 

DELT 

DEN 

DENI 

DEN2 

DENS 

DFLO^ 

DFIN 


J-S, used in computing interface location 

value of integer M converted to floating point 

length of section, ft 

volume expansivity of wax, R"^ 

J-2, used in computing interface location 

thermal capacitance of node (I,J), Btu/°F 

constant pressure specific heat of wax, Btu/lbm°F 

constant pressure specific heat of wax, Btu/lbm°F 

constant pressure specific heat of bottom plate, Btu/lbm°F 

constant pressure specific heat of fin, Btu/lbm°F 

constant pressure specific heat of top plate, Btu/lbm°F 

average height of liquid based on amount melted, in. 

interfacial velocity for nodes adjacent to centerline (I = MM), in./h 

absolute value of temperature difference between bottom plate and 
interface, °F 

wax density, Ibm/ft^ 

bottom plate density, Ibm/ft 

fin density, Ibm/ft^ 

top plate density, Ibm/ft^ 

approximate interfacial location for nodes adjacent to fin (I = 2) based on 
amount melted being equal to F2JM, in. 

interfacial position for nodes adjacent to fin (1= 2) based on any amount 
being melted, in. 
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DIF 

DMID 

DMO 

DT 

EPS 

ERROR 

F(U) 

FAC 

F2JM 

G 

H 

HBOT^ 

HMELT 

HT0P3 

HTR 

I 

J 

JOE 

KCHK 


temperature difference used in comparing new and old temperatures 
during iteration when solving steady state equations for unspecified fin 
temperatures, °F. 

interfacial position for node adjacent to centerline, in. 

DMID evaluated at previous time, in. 
time increment, h. 

arbitrarily set small number used as a comparator 

percent error in computed energy balance based on transfer rates, percent. 

mass fraction of node (I,J) which has undergone phase change since start 
of process 

time ratio used in linearly interpolating specified fin temperatures at a 
particular time in terms of bracketed data values. 

fraction of S which corresponds to 1/32 in. (arbitrary) 

acceleration of gravity, ft/h^ 

PCM section height (see Fig. 40), ft. 

heat transfer coefficient between external fluid and bottom plate, 
Btu/h-ft^-°F 

heat of fusion, Btu/lbm 

heat transfer coefficient between external fluid and top plate, Btu/h-ft^-°F 
heat of transition, Btu/lbm. 

integer designation of vertical column in which a mode is located (see 
Fig. 40) 

integer designation of horizontal row in which a mode is located (see 
Fig. 40) 

counter used in refining the heat transfer computation before progressing 
in time 

integer used to control printing of results at desired times (see definition 
of KCOUNT) 
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KCOUNT 

M 

MCOUNT 

MFIN 

MM 

N 

ND 

NDP 

NI 

NJ 

NN 

PR 

Q(I,J) 

QBT 

QBW 

QFTR 

QMELT 

QRAT(U) 

QS(I,J) 

QSIN 

QSW 


integer counter used to print our results at times when KCOUNT = KCHK 

number of wax nodes in a horizontal row 

counter used in determining unspecified fin temperatures 

maximum value of MCOUNT which when exceeded causes program to 
stop 

M+1 (see Fig. 40) 

number of wax nodes in a vertical column 

number of data points for measured fin temperatures 

ND-1 

N+4 

N+3 

N+2 

Prandtl number 

unnecessary variable -- replaced where needed by QS(I,J) 

instantaneous heat transfer rate through bottom, Btu/h 

instantaneous heat transfer rate through bottom to wax only, Btu/h 

instantaneous rate of heat transfer to fin, Btu/h 

energy which accounts for amount of wax melted at any time, Btu 

the instantaneous net rate of heat transfer to node I,J, Btu/h 

the energy stored by node 1,J above TREF for wax and above 0 for metal 
nodes, Btu 

the energy stored by node I,J above TREF corresponding to initial 
temperature throughout network, Btu 

instantaneous heat transfer rate from fin to wax, Btu/h 
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TIME 

TIN 

TK 

TKL 

TKR 

TMELT 

TREE 

TTR 

T1(I,J) 

T2(I,J) 

T3(1,J) 

TAU2 

TKI 

TK2 

TK3 

TIM(L) 

TMl(L) 

TM2(L) 

TM3(L) 

TM4(L) 

VIS 


instantaneous value of time, h. 
initial temperature of all nodes, °F 

wax thermal conductivity (artificially allowed to vary in liquid to account 
for convection), Btu/h-ft-°F 

thermal conductivity of liquid, Btu/h-ft-°F 

ratio of effective thermal conductivity to thermal conductivity 

fusion temperature of wax, °F 

arbitrary reference temperature (should be less than TIN), °F 
transition temperature, °F 
temperature of node I,J at time t, °F 
temperature of node I,J at time t+At, °F 

temperature of node 1,J in fin at beginning of each iterative step used in 
finding steady state solution, °F 

arbitrarily defined time value used in print-out control, h. 
thermal conductivity of bottom plate, Btu/h-ft-°F 
thermal conductivity of fin, Btu/h-ft-°F 
thermal conductivity of top plate, Btu/h-ft-°F 

time value corresponding to input data of measured fin temperatures, h 

measured bottom plate temperature (input data), °F 

first measured fin temperature (input data), °F 

second measured fin temperature (input data), °F 

third measured fin temperature (input data), °F 

viscosity of liquid, Ibm/h-ft 
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VI volume of wax melted at time t, ft® 

V2 volume of wax melted at time t+At, ft® 

W width of wax cell, ft 

1. This assumes that some finite thickness must have melted before it could be 
detectable on the film. The number DFIN is the height corresponding to a node with 
any amount melted. 

2. This was included to be general but has not been used to date as bottom plate 
temperatures were specified as input data. 

3. This has been included but set at a small value to essentially correspond to the top 
being insulated. 
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D.3. DISCUSSION OF COMPUTER PROGRAM FOR MELTING 


In the following discussion, references are made to line numbers corresponding to 
those designated on the copy of the program. 


LINES 

1-4 


5-23 

24-54 

55-63 

64-82 

83-95 

96-108 

109-112 

113-144 


145-150 

151-172 

173-174 


DISCUSSION 

required dimension statements for subscripted variables; values should be 
(MM, NI) for all double subscripted variables except T3 for which they 
should be (1, NJ); values should be ND for single subscripted variables. 
NOTE: Q(MM,NI) is superfluous and can be omitted with lines 136-140. 

specification of computational parameters 

specification of physical properties 

specification of geometry parameters 

computation of all vertical thermal resistance values RV (1,J) 

computation of all horizontal thermal resistance values RH(l,J);note that line 
95 is a special definition which amounts to bypassing node (l,N+3) which 
was done to overcome stability criterion required by this small corner node 

computation of all nodal capacitance values C(1,J) 

computation of energy stored by a wax node relative to TREE for the start 
and end of phase transition and the start and end of fusion, respectively 

initialization of pertinent quantities; the initial value of the stored energy 
depends on the relationship of the initial temperature to the reference 
temperature TREE; note that lines 136 140 are superfluous and can be 
omitted; some initialization of certain parameters is done in the DO loop 
between lines 160 and 172 which are set primarily to avoid random printout 
and are not essential to the computation done in the heart of the program 

input data values foi measured fin and bottom plate temperatures and 
corresponding time values 

printout of initial values for checking purposes and printout of certain 
computed quantities for informational purposes 

initially defines T3(1,J) for all fin nodes and sets these equal to the initial 
temperatures T 1 ( 1 ,J ) 
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175 beginning of main computation scheme 

176 time is stepped forward by At 

177-188 loops which compute and store the net rate of heat transfer to nodes (1,J) 
based on old temperatures T1(I,J); the rate of heat transfer to node (I,J)'is 
given by 

ntl Jt = T(I-1,J)-T(U) ^ T(I-H,J)-T(I,J) ^ T(1,J-1)- T(I,J) ^ T(I,J+1)-T(U) 
RH(I,J) RH(I,J) RV(I,J) RV(I,J+1) 

this expression must be modified accordingly for nodes near a boundary 
which are not surrounded by four neighbors 

189-191 computation of total energy stored by node (I.J) since the start which is 
given by 


^stored 


tUT.e 


192-193 computation of new top plate temperatures T2(I,N+3) from the expression 


T2(I,J) = T1(I,J) + 

C(I,J) 

194 sets the corner plate node (l,N+3) temperature equal to that of the second 

node (2,N+3); omission of the corner node in the computation scheme was 
done to avoid stability problems due to its small size 

195-208 computation of new wax node temperatures from the energy stored by the 
nodes and their capacity and/or phase change enthalpy values; when the 
stored energy lies between Q1 and Q2 the new temperature is forced to be 
the transition temperature and when it lies between Q3 and Q4 the new 
temperature; is forced to be the fusion temperature; also the fraction of the 
node which has undergone phase change is calculated from the relationship 
of the stored energy to Q3 and Q4; note that certain designated cards need 
to be changed when running the program for freezing as contrasted to 
melting 


209—210 beginning of determination of fin temperatures 
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211 - 
213 - 

221 - 

223 

224 


238 

239 

240 

243 

255 

258 


212 all new fin temperatures are set to T3(1,J) which simply represents an 

assumed value always corresponding to the previously computed value except 
at the very beginning at which time it is set as the initial temperature 

220 interpolation scheme which assigns new temperatures to the three nodes on 

the fin and one on the bottom plate corresponding to positions where 

temperature measurements were made; the new temperatures are linearly 
interpolated from the input data 

-222 assigns all nodes along the bottom plate the same value of new temperature 

beginning of iteration process to determine unspecified fin temperatures from 
steady state equation; fin nodes were not treated as transient cases due to 
their extremely small capacitances that would impose a severe stability 
criterion 

-237 iteration process used to determine unspecified fin temperatures; in each 
iteration, T2(1,J) is computed from steady state equations and then 
compared with T3(1,J) which corresponds to the calculated temperatures 
during the previous iterative step; the interation is continued until the 
differences between computed fin temperatures and their corresponding 
values for the previous iterative step are all acceptably small; should the 
iteration exceed MFIN counts the program is directed to stop 

when counter JOE equals 2, the new temperatures at all nodes are 
considered to be the solution at the particular time and the program 
advances to line 259 

increase of counter JOE to 2 

_242 the energy stored at each node (I,J) is reset back to its original value, this is 
to allow for an improved computation of the net heat transfer rate to each 
node to be made and then a recomputation of the stored energy and the 
corresponding new temperatures 

-254 computation of net heat transfer rate to each node using new temperatures 
T2 

-257 calculation of net heat transfer rate to each node as the average of that 
based on old temperatures T 1 and new temperatures T2 

return to line 189 which consists of redirecting the computation through 
that of computing improved new temperatures, energy storage values, and 
fractional melted values by using the improved (averaged) heat transfer rate 
(Lines 255-257); this corrective technique is only employed once 
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259 

279 

283 

284 
285- 

287 

288 
289- 
291- 

294 

295- 

316 

317 

318 
319- 
321- 

333- 

336 

338- 


278 computation of various heat transfer quantities from the new temperatures 
obtained at time t+At 

-282 computation of volume of melted wax 

computation of average liquid depth from the volume melted 

computation of energy required to melt the wax which has melted 

-286 resetting of VI and V2 for next time step 

determination of absolute value of temperature difference between bottom 
plate and the fusion temperature 

computation of the Prandtl number 

-290 computation of the Rayleigh number 

-293 determination of effective liquid conductivity due to convection by using 
correlations of O’Toole and Silveston 

computation of ratio of effective liquid conductivity to actual value 

-315 recomputation of thermal resistances in the liquid by using the effective 
thermal conductivity rather than the actual value 

counter check which controls printing out of desired results as well as 
computation of interfacial position 

counter advance 

by-pass of printing results except when line 316 is executed 
-320 write statement for printing results 

-332 computation of interface position for columns next to fin and next to 
centerline 

■335 write statement for printing results 

337 calculation of interfacial velocity and renaming interfacial position to provide 
for determining its change at the next time step 

351 write statements and corresponding formats for printing results 
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352- 

359- 

361 

363 

364 

366 

367 


-358 reinitialization of pertinent quantities for next time step 

361 setting new temperatures for current time step to be old temperatures for 
the next time step 

reinitialization of counter 

comparison of time to upper limit value which when exceeded results in 
stopping the program 

-365 printing out of value of counter used in fin temperature iteration 
STOP 
END 
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Figure D-1 . Skelton flow chart. 
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D.4. COMPUTER PROGRAM FOR SOLIDIFICATION 

TEST 230-15 FREEZE TEST 3/4 INCH CELL 

DIMENSION RVI5#31)# RH(Si31)^ CI5#31)/ 08(5^311/ Tl(5 
1 « 31 ) 


dimension T2(5>31)* F(5/3lJi QRAT(5/31)# T3 ( 1/ 31 ) / QRA2 
i ( 5# 31 ) 

dimension TMK25)# TM2(25li TM3(25)/ TM4I25)/ TIM(25) 
dimension VR(5*31)/ HR(5<3t) 

c computational parameters 

N-27 


M-4 

AM-M 


MM-M+1 

NN-N+2 

NJ-N+3 

NI-N+4 

ND-23 

NDP-22 

TAU-li5 

KCOUNT-1 

MCOJNT-l 

MFlN-500 

EPS-l»E-06 

KCH<-1000 

JOE-1 

DT-0»00005 

TAU2-I 10«»DT )+{DT/3. ) 

- AF-4Qi 

UN-1»J-EPS 

C PHYSICAL PROPERTIES 
TAM3-80* 

TIN-lOOt 
HTOP-1 .£-08 
HBOT-5. 

Q-02«174«3600i*3600« ) 
C WAX 

0En-47»2 
TK-0»087 
CP-0.5 
TTr-73 * 04 
HTR-22* 108 
TMELT-89.8 
HMELT-73.357 
TREF-30. 

BETA-0.00045 

VlS-14.3 

CL-CP 

tkl-tk 

C BOTTOM PLATE DENOTED BY I 
DENl-171. 

T<l-93. 

CPl-0.22 
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c TIN DENOTED BY 2 
DEN 2-171 • 

T< 2 - 93 » 

CP 2 » 0»22 

C TOP PLATE DENOTED BY 3 
DEN 3 » 72»5 
TK 3 » 0*09 
CP 3 - 0 O 3 

c geometry parameters 

H- 0 » 75 / 12*0 

H- 2 » 625 / 12*0 


C 


C 


c 


B«5»0/12»0 
S-W/ ( 2 * »AM ) 

Sl»0»032/12*0 

S2«0»008/12»0 

S3«0»25/12»0 

F2JM-1 'O/l 32'0^S412*0 ) 

vertical resistances 

RVJli2)»(Sl/(TKl»S2*0))+(2«/l H0OT*S2*B ) ) 

DO 10 I-2JMM , 

10 RV(I/2I"(S1/12»»T<1*S»B) )+(1»/(HB0T*S*B) > 
RV{ l/3)-(S/( TK2*S2*B) )+<Sl/(TKl*S2*B) ) 

DO 20 I»2iMh 

20 RVIIi3|-(l»/(2*#T<*B>) + ( Si/ ( 2* ¥Tk 1*S*B ) ) 
DO 30 j»AiNN 

30 RV(l/J)-( (2*»S)/(T<2»S2»S) ) 

DO HO NN 

DO HO 1-2/MM 
HO RV( I#J)»( 1*/(TK»B) ) 

RV{liN + 3)-(S3/lTK3*S2¥0J )>(S/< TK2*82*B J ) 


( 1 »/(HT 0 P*S* 8 ) ) 


DO 50 1-2/MM 

RVM/N + H) - (S3/(2«*TK3*S*6M + 

50 RV(I/N + 3)-(S3/(2*¥TK3*S*B)) + ( I •/< 2**TK*B > ) 

RV(1/N+H)-(S3/(T<3*S2*B) )+(2«/IHT0P*S2*B) ) 
RV(2/N+H)-(S3/(2»»Tk3¥0*<S+(S2/2« ) ) ) ) ♦ { 1 • / < HTOP*B* < S+ 

1 J S2/2» 1 ) ! ) 

horizontal resistances 

DO 70 1-2/MM 

RHtI/2) - t S/ ( T<1 »Sl*B I > 

70 RH( I/N+3»-(S/(TKJ¥S3*3) ) 


DO 80 J-3/NN 
00 80 

80 RH ( 1/ J ) - ( 1 • / ( TK*B ) 1 

DO 90 J-3/NN , ^ ^ 

90 RH(2/J)« (S2/(2»i‘T<2‘i‘S'^B) ) ♦ ( I • / ( 2 • *T<»B ) ) 

RHl2/2)-( (S2 + S)/(2.»TKl*Sl*B) ) 

RH(2/N + 3)-( (S2 + S)/( 2« »TK3*S3*B ) ) 

RH(2/N + 3)-RH(2/N + 3)+RVU/N + 3) 

nodal capacitances 

C( 1/2)- (( denies 1»S2*B*CP1)/2* ) 

DO 110 1-2/MM 
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c 


110 C ( I* 2 ) ■ ( DEM^Sl ¥S»B»CP1 ) 

DO 120 J« 3 iNN 

120 C(l/J)«( (D£N 2 *S 2 »S<^B^^CP 2 )/ 2 *) 

C ( 1 < N + 3 ) ■ ( { DEN 3 <‘S 2 ¥S 3 ¥B»CP 3 ) /£• ) 

DO 130 I« 2 /MM 

130 C { I # N+ 3 ) ■ ( DEN 3 »S 3 »S*B»CP 3 ) 

C( 2 #N + 3 )«(S 3 afB*(S+(S 2 / 2 . ) )¥DEN 3 *CP 3 ) 
DO 140 J« 3 iNN 
DO 140 I- 2 /MM 


140 C ( I # J ) ■ ( 0£N¥ ( S<‘¥2 ) »B*CP ) 

Q1 ■ (C(2i3)»( TTR-TREFI ) 

Q 2 »ai + ( ( 0 EN¥(S 4 '¥ 2 )i^B‘i‘HTR)) 

03 - 02 + ( (C(2/3) )¥(TMELT-TTR) J 

Q4»Q3+( (DEN¥(S»¥2) ¥&¥HMELT ) ) 

Initialization of pertinent quantities 
time»o»o 

QKAX»0*0 
QBn > 0 • 0 
OSWaO • O 
QTw»0 • 0 
OTOP-0.0 
QFTR»0»0 
QTTRaO • 0 
V1«0#0 


V2-0»0 
DMO»0»0 
DO 150 J«2iNj 
DO 150 I-l/hM 
150 Tl( I/J)-TIN 


IF(TlN«LT*TTR)QSlN«(C(2i 3 )♦ « TIn-TREF ) ) 
lijTRM * 08IN«Q2+(C(2>3)*| TIN 


IF(TlN.eT.TMELT)QSIN-Q 4 +nC( 2 > 3 ) )»(TIN-TMELT) ) 

DO 160 I«l>MIi 
160 Tl(IiN+4)aTAMB 
DO 170 J"3 >Nn 
DO 170 I*2rfhh 
170 OS(I>J)-qSIN 
DO 199 J-2,M 
DO 199 I-l>hM 
F ( I / J ) aO • 0 
199 QRaT ( li J ) aO.O 

READ<5/11) (Tril(I)i Ial,NO) 

READI5/11) (TM2(I)i lal^NO) 

READ(5>11) (TM3(I)i Ial,NO) 

REaD( 5/11) (TN4(I), Ial,NO) 

REA0(5ill) (Tlf-KI)/ lal^NO) 

11 FORMAT! 8F10»0 ) 

DO 910 KTC-li ND 
910 TlM(KtC)aTlM(l<TC)/ 3600 » 

WRITE(6^22) TImE,w/M>N/M 
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28 FORMAT ( lXi5MTIiit«iE15»8*10X/?Mfi(«iE15»8*10Xi2HH*AEl5*8 

1 # 1 OXi 2HN»/ I E# 1 5X I 2 HM»j I 2 ) 

WRITE{6>33)Ql>Q2^03iQ4 ^ „ 

33 F0RMATIlXj3HQl*#E15*8il0Xi 3H02»#El&»8# 10X^3HQ3*aE15«8 
l/10X/3MQ4»iEl5»8) 

WRITEI6/44) 

44 F0RMAT(2XMH1 J/ 5Xi ISHVERTICAl RESISTANCE#6X 

I# 21HHORIZONTal RES1STANC£#6X# 17HNOOAL C APaC I T AMCE* 6X 
2/11hTEmPERATuRE> 5X#7HQS(I>J)I 
DO 200 j»EiNl 
DO 200 

IF ( J»EU»N+4 ) GO TO 1 
GO TO 2 

1 QS< Ii Jl-OtO 
RH( I# J)-l *E08 
C ( li J ) -O^O 

2 IF(IiEQ»l) RH ( I i J ) "I • E08 
IF<J*EQ*2«0R*J‘EQ*N+3) QS(I^J)»OtO 
IF(I.Eu.l) QSlIiJ)»0*0 

wRITE(6#55)IiJiRv(IiJ)iWH(IiJ>#C(I#J)#TlU#J)#Q9(I#J) 

55 format t lX> 12/ IX/ I2/5X* E15»8i ICX* E15«8# lOXiE15.B^6X 
1/E15*8/4X/El5*ai 

200 continue 

DO 400 J»3/NJ 

DO 400 I«2/MM 

VR(I/J)»RV1I/J) 

aOO HR( 1/ J)»RH1 1/ J) 

DO 889 J-3/NN 

aS9 T3( 1/ J)-T1 ( 1/ J) ^ 

c computation section computation section computation 
c SECTION 

3 TIME-TIMEi-DT 

do 210 J-3/NJ . ' . , 

QRAT(l/J)-i(Tl(2/J)-Tl(l/JJ 1/RH(2/J) I ♦( ( T 1 ( 1 / J-1 ) -T I ( 1 
1/J) )/RV(l/Jl (T1(1/J + 1)*T1(1/J) J/RVll/J + l) ) 
QRAT(MM/J)“t (T11M/J)“T11MM/J) )/RH< MM/ J) ) + l (TKMM/J-l) 
1-T1(MM/J))/RV(MM/J))-M(T1(MM/J4-1)-TI( mm/ j ) ) / RV ( MM/ J 
2 + 1 ) ) 

210 continue 

IF( M.EQ* 1 )G0 TO 212 
DO 211 J-3/NJ 

DO 211 1-2/M ^ ,, ,, 

QRAT( 1/ J)-( (T1 ( 1-1/ J)-T1( I/Jn/RH(I/J) ) + MTU 1 + 1/ J)-Tl 

1(I/J) )/RH( I + l/ J) ) + l (Tl(I/J“l)-Tl(I/J) )/RV(I/j) )♦( (Tl(I 
2/ j+1 ) -Ti ( 1/ J ) I /RV ( 1/ j+1 n 

211 continue 

212 DO 220 J-3/NJ 
DO 220 I-l/MM 

220 OS<I/J)-QS(I/Jl+(QRAT(I/J>+DT) 

DO 240 I-2/MM 

240 T2tI/N + 3)-TKI/N + 3) + ( (QRATC 1/ N + 3 ) +DT I /C ( I/N + 3) ) 
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non on 


T2( l>N+3)«T£(2>N+3) 

DO 250 J«3iNN 
DO 250 I«2>NM 

IFOS(IiJ).LT*01) T2(I>J) ■TREF+(0S( li J)/C ( 1/ J) ) 
IF{asn/J).QE»Ul.AND*OS(IjJ)»LE» 02 ) T2(I>J)»TTR 
IF(3S(IiJ)»GT*Q2«AND*QS(I#J)*LT«03) T2|I#J!«TTR+((QS(I 
l/J)*Q2)/C(IiJ) ) 

IF(QS(liJ)«GE»Q3»A\D*QS(I/J)tLE»Q4) T2(IiJ)«TlELT 
C WHEN aOiNQ FROn melt TO FREEZE OR ViCE-VERSA CHANGE THE 
C FOLLOWING CA 

IF(3S( I/J) •GT»03«AND*QS( I#J) »LT»Q4) F( I/JJ*I04*QS( I 
l/JJJ/iDEN^lS'M'El^B^HMELT) 

IF(3S( 1/J) .GE.J4) F(I,J)«0«0 

IF(QS(I>J).QT*g4) T2( li J >»TMELT+( ( QS ( I # J ) -Q4 > /C ( I < J ) ) 
WHEN GOING FROM MElT TO FREEZE OR VICE-VERSA CHANGE THE 
following CA 

IFOSI 1/ J) .LE»Q3) F(I/JJ«1.0 
250 continue 

specification AND/OR determination OF fin temperatures 
The following do loop assumes fin temperatures 
FOR iteration 
00 255 J-3iNN 

255 T2(l/J)-T3( 1/J) 

00 260 L-1#NDP 

IF t TIME »GE • TIM ( L )• and. TIME 'LE* T iM ( L+i ) ) GO TO 4 

260 continue 

4 FAC»(TImE-TIM(L) )/(TIM(L + 1 )-TlMlL) ) 

T2(l/2)»TMl(L)+( (TM1(L+1)-TM1 (l) ) ♦FAC ) 

T2(1#8)«TM2(L)+( (TM2IL+1J -TM2 I l ) ) *FAC ) 
T2(lii5)»TM3(L)+((TM3(L+l>«TM3(L) )»FAC) 
T2(1#21)«TM4(L)+((TM4(L+1I-TM4(L))»FAC) 

DO 270 I»1/MM 
270 T2( I#2)-T2( 1^2) 

C UNSPECIFIED FIn TEMPERATURES DETERMINED BY STEADY STATE 

c equations 

256 MCOJNT-MCOUNT+1 
DO 280 J-3jNN 
T3(1/J)-T2(1>J) 

IFl J.EQ.8»0R.J.EQ.15) go TO 5 
1F(j.EQ.21) go to 5 

T2( 1, j).( ( T2( J-1 )/RV( li J) )+( T2(2# J)/RH«2#U) )+<T2l 1/ J 
l + l)/RV(l/J + l) ) )/( (l./RVdiJ) )♦( I ./RM( 2/U > ) ♦ ( I »/RV ( 1# J 
241 ) ) ) 

5 continue 

280 CONTINUE 

IF( HC0UNT.QT.MFIN ) GO TO 3 
DO 281 J-3iNN 
DIF»T2( li J )-T3( 1/ J) 

IF( ABS( DIF ) .GT.EPS ) GO TO 256 

281 continue 

IF(jOE.EQ.2) go to 285 
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• I )-T2< I 


jOE-JQE+1 
DO 88d J«3iNJ 
DO 888 I-l#hPl 

838 QS( I#J)“QS( 3RAT ( IiJ)*OT) 

QRA2UiJ)«( 

IjJ) >/RV(1iJ) )+( )/RV( l/J+l) I 

QRA2(MMiJ)-( (T2(M*J)-T2{MMiJ) 

l-T2(M^#J) )/RV(rlM/ J) ) + ( (T2(MM# J + l I-T2JM1# J) )/RV(rtf1# J 
2*i ) ) 

282 CONTIMUE 

IF( ?1*Ea» 1 IGO TO 283 
DO 283 J»3iNJ 

Q°A2a^ j”-nT2( I-i# JI*T2(Ii J) »/RH( 1/ j) ) + MT2( 1 + 1/ J)-T2 

i?ujM'RiilU:jnH(Taa.j.i.-T2(t.jn/Rvn.jM*nT8(i 

2iJ + l)-T2( I/J) )/RV( I^J + l) ) 


283 CONTIMUE 

DO 284 J«3>NJ 
DO 284 I-lihM 

284 QRATI I# J)-1QRAT( li J)+QRA2( liJ) )/2»0 

QO TO 212 

285 MCOUNT-1 

DO 401 J»3<NJ 

DO 401 I«2iWh rn iLft a~ 

IF ("F ( J ) »GT »EPS » AMD»F( I< J) *LT«UN ) GO TO 402 

IF ( F ( I< J ) * GT « UN ) 30 TO 403 

OO TO 401 — 

403 RV( I* J)»VR< 1/ Jl 

RH(IiJ)-HR(IiJ) 

QO TO 401 

»02 RV ( I i J ) »VR < J ) /AF 

RH( I> J)«HR( 1/ J)/AF 

RVt li Jt-1 )«VR( li J-t-l )/AF 

RHt l4l/J)«riR{ I*l> J)/AF 

401 CONTIMUE 

DO 286 I»1/MM 

286 QRAT( I^2)«( (T2(I>21-T1(I<2) l♦ClIJ2)/DT) 


287 

290 

300 


DO 287 MM 

0S{ I>2)-QS( I/2)+(C( I#2)a<T 2( I* 2)-Tl( 1/2) » ) 

DO 290 1-2/MM 

0BW-(T2(I/2)“T2{I/3) ) /RV ( 1/3 )+3BW 
qTW-( T2( I/N + 3) -T2( 1/ M + 2 J >/RV ( I»N + 3)+QTw 
qBT-Q8w+{ (T2( 1/2)-T2(1/3) )/RV( 1/3) ) 

DO 300 J-3/MM 

QSvy«QSW+( (T2( 1/J)-T2(2/J) )/RH(2/J) ) 


DO 310 I-l/MM 

310 qTOP-( (T 2(I/Nt3)-TAMB)/RV<I/N+4) )+OTOP 
DO 311 J-3/MN 

311 QFTR-QFTR+QRaT ( 1/ J ) 


DO 312 I-l/MM 
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312 OTTR-QTTR+QRaT{ liNJ) 

OWAX- r QBW + OTW + QSW ) ¥DT+IQWAX ) 

RATIO. qsw/OBW 

ERROR.( (aBT-(QSW + OSH4QTK + QTOPi-3FTR>QTTR) U100«0)/QBT 
DO 930 J-3/NN 
DO 930 I.2/MM 
V2«V2+(F(I^J)»(S¥¥2)¥B) 

930 continue 

DAVQ.( (2»0¥V2)/(wi‘BI l¥l2»0 

OMELT»( (V2-Vi ) ¥DEN¥HhELT)/DT 

QFREZ-UhELT 

V1-V2 

V2«0*0 

IF(<C0UNT»EQ*KCH< ) 00 TO 6 

KC0UN7»KC0UNT+1 

QO TO 7 

6 WRITE(6/66) TIME# OWAX# ERROR 

66 FORMAT! 1 X# 5HT I ME.# E 1 5 • 8# lOX# 5H3WAX-/ El 5 • 8# lOX# 6HERR0R- 
1#E15*8) 

DO 313 J-3#NN 

AJ-J-3 

BJ.J.2 

IF(F(2#J).gT»F 2JMI DFL0«(BJ¥S*12*0 J 
DFIN«( AJ»S¥l2»0) 

IF(F(2# J) .LT'EPS) GO TO 314 

313 CONTINUE 

314 DO 315 J-3#NN 
AJ-J-3 

IF(F(nM#J) •LT»1»0) DMIO«( (Aj¥S)-f(Fjf1M#j)¥S) )¥l2»0 
IF(F(!1M#J) *LT«liO) GO TO 316 

315 CONTINUE 

316 WRITE<6#67) DFIN# DFLO# OAVQ# OMID 

67 FORMAT! 1 X# 5HDF I N«# E 1 5 • 8# lOX# 5h0FL0-i E15 • 8/ lOX# 5M0 A VQ« 
1#E15»8# 10X#5H0MI0«/E15»8) 

DDOT"!DMID-DMO )/DT 
DMO-DMID 

WRITE!6#69) QFREZ# ODOT 

69 FORMAT! 1X#6HQFREZ-#E15«8#10X#5HDD0T-#E15«8) 

WRITE!6/77) QBW# 3SW# QTWi QBT# RATIO 
77 FORMAT! 1X#4 HQBh.#£15»8#3X#4HQSh(«#E15»8#3X/4HQTW"/E15*8 
1#3X#4HQBT«#E15«8#3X#6HRATI0«/E15*8 ) 

WRITE!6#88) 

88 F0RMAT!2x#1HI#2X#1HJ# 1 OX# 1 IMTEMPERaTURE# 10X#15HFRACTI0 
IN MELTED# lOXi 1 IHENERQY RATE# IOX# 13HENERQY STORED) 

DO 320 J»2#NJ 
DO 320 I.1#MM 

wRITE(6#99) I# J#T2! I# J)#FI I# J)#3RATI I# J)#QS( I# J) 

99 FORMAT ( IX# 1 2# 1 X# 1 2# 8X# E 15 *8# 8X# E15 • 8# 8x# 1 1 5 • 8# 8X 
1#E15*8) 

320 continue 
KCOUNT- l 
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7 QBW»0»0 
QSVt*0*0 
QTw«0»0 
qTOP»0»0 
OTTR»0*0 
qFTR-0.0 

DO 330 J*2#NJ 
DO 330 

330 Tl ( li J)«T2( 1/ J> 

J0E«1 

IF ( TIME »LT »TAU) QO TO 3 

8 WRITE(6#111) MCOUNT 
111 FORHATt IXi 13) 

STOP 

END 
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